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INTRODUCTION 


The  Naval  Facilitiea  Engineering  Command  seeks  to  expand  the  tech- 
nology  base  upon  which  future  shore  facilities  are  founded.  An  important 
area  in  which  a  new  design  approach  is  being  developed  is  in  the  nanage* 
ment  and  understanding  of  those  characteristics  of  earthquake  and  explo¬ 
sively  generated  foundation  motions  that  damage  or  destroy  Naval  equipment. 
This  report  is  a  part  of  the  design  method  development.  Specifically, 
it  documents  and  offers  a  new  method  of  shock  and  response  spectrum 
computation.  The  new  method  is  more  accurate  than  any  known  to  be  in 
use,  it  runs  on  a  computer  with  less  computer  time  than  other  procedures 
with  similar  sophistication,  and  the  theory  is  more  easily  derivable  and 
understood.  The  shock  spectrum  is  the  concept  used  in  the  new  design 
method  for  quantifying  the  destructive  capacity  of  explosively  and 
earthquake  generated  violent  equipment  foundation  motions. 

The  report  begins  with  an  overview  of  the  shock  and  response  spec¬ 
trum  concept,  and  then  describes  the  computations  that  are  required  to 
transform  time  histories  into  shock  spectra.  It  then  discusses  current 
computation  methods  and  compares  them  to  the  new  computation  method. 

The  new  method  of  computing  the  "during"  values  of  the  spectrum  is 
presented  here  for  the  first  time;  the  theoretical  detail  is  given 
separately  in  Appendix  A.  The  procedure  for  the  computation  of  the 
residual  spectrum  values  has  been  presented  previously  in  an  interim 
report  (Ref  1)  but  is  repeated  here  for  completeness.  A  FORTRAN  IV 


listing  of  the  programs  is  given  in  Appendix  B,  as  it  is  no 
on  a  time-sharing  computer. 
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BACKGROUND 


The  term  “shock  spectrum"  has  been  in  use  for  about  40  years,  and 
so  one  might  reasonably  assume  the  technology  to  be  complete.  Respected 
firms  now  sell  preprogramed  minicomputers  that  calculate  and  plot  a 
shock  spectrum  at  the  touch  of  a  button.  Thus,  it  is  felt  that  some 
explanation  for  yet  another  report  on  the  subject  must  be  offered. 

A  precise  definition  of  the  term  “shock  spectrum"  has  not  yet  been 
accepted.  The  term  is  used  by  DOD  and  itB  contractors,  while  the  earth¬ 
quake  community  synonomously  use  the  term  “response  spectrum."  The  term 
is  used  mostly  to  describe  a  short-duration  violent  motion,  but  it  also 
is  used  to  describe  a  force  transient.  For  the  case  where  it  is  used  to 
describe  a  motion,  it  is  invariably  "defined"  to  be  the  peak  response  of 
single-degree-of-freedom  systems  to  that  motion  plotted  versus  the 
natural  frequency  of  the  single-degree-of-freedom  system.  Within  that 
definition  there  are  at  least  108*  different  plots,  even  all  with  the 

same  damping,  which  one  could  call  the  shock  spectrum  of  that  motion. 

This  state  breeds  confusion,  and  prevents  one  from  acquiring  experience 
in  the  appearance  of  severe  shock  spectra. 

Thus,  a  precise  definition  is  required  for  use  in  any  design  proce¬ 
dure,  and  one  was  given  along  with  reasons  for  it  in  the  preliminary 

design  method  (Ref  2).  Specifically,  the  shock  spectrum  is  taken  to  be 

a  plot  of  the  peak  relative  displacement  of  a  single-degree-of-freedom 
system  exposed  to  the  motion  being  analyzed,  as  a  function  of  undamped 
natural  frequency,  and  plotted  on  four-coordinate  paper.  The  shock 
spectrum,  therefore,  is  treated  as  a  precise  technical  property  of  a 
transient  motion.  A  property  is  defined  by  how  it  is  measured  and, 
particularly  in  the  case  of  shock  spectra,  by  how  it  is  actually  computed. 
This  report  in  that  sense  defines  the  shock  spectrum. 


♦Plot  it  log,  linear,  or  semi-log,  or  four-coordinate;  plot  absolute  or 
relative;  plot  peak  acceleration,  velocity,  or  displacement;  give  plus 
and  minus  values,  or  the  overall  value;  plot  the  during,  residual,  or 
the  total  peak.  Four-coordinate  is  the  same  as  log  but  the  paper  has 
more  lines. 
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Although  defined  and  specified  for  purposes  of  this  work,  interest 
and  additional  study  of  the  shock  spectrum  will  continue.  Analysis  of 
dynamic  data  in  terms  of  its  effect  on  single-degree-of-freedom  systems 
is  a  basically  fundamental  approach  that  has  not  yet  been  fully  exploited. 
The  references  for  this  work  are  current;  work  on  explaining  the  calcu¬ 
lation  and  improving  accuracy  continues.  Appendix  A  of  this  report  is  a 
new  development.  The  comments  at  the  end  of  Appendix  A  indicate  that 
even  this  effort  can  be  continued.  However,  the  current  state  of  shock 
spectrum  understanding  is  adequate  for  a  useful  design  method.  This  is 
explained  in  the  preliminary  design  manual,  Reference  2.  For  those  who 
have  to  compute  their  own  spectra,  or  use  preprogramed  machines  and 
desire  an  appreciation  for  the  computations,  this  report  should  help. 

The  shock  spectrum  is  the  key  to  an  upgraded  design  method  for 
installation  of  shock- resistant  equipment.  Pursuance  of  the  concept  can 
make  facilities  much  safer  for  Navy  personnel.  For  example,  a  common 
problem  from  earthquakes  is  falling  fluorescent  light  fixtures.  Attention 
to  the  shock  spectrum  can  cause  one  to  use  a  fixture,  not  necessarily 
more  expensive,  that  is  not  sensitive  to  motions  with  earthquake  fre¬ 
quencies.  Indeed  all  appurtenances  to  the  structure  and  all  installed 
equipment  have  estimatable  natural  frequencies.  The  shock  spectrum 
gives  the  maximum  values  of  oscillatory  motion  to  expect  from  equipment 
with  various  natural  frequencies.  Dynamic  design  is  considerably  sim¬ 
plified  when  one  can  predict  the  peak  acceleration  and  displacement. 

Thus,  the  mere  understanding  of,  and  attention  to,  the  information  in 
the  predicted  shock  spectrum  permits  an  evaluation  of  installed  equipment, 
especially  as  regards  the  equipment  breaking  free  and  injuring  personnel. 


DESCRIPTION  OF  THE  SHOCK  AND  RESPONSE  SPECTRUM 

The  explosively  generated  mechanical  shock  motion  or  the  earthquake¬ 
generated  motion  of  equipment  foundations  is  often  recorded  as  a  signal 
from  an  accelerometer  on  magnetic  tape.  These  records  are  converted  or 
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digitized  into  a  list  of  closely  spaced  numerical  values  of  the  acceler¬ 
ation  with  known  equal  time  intervals  separating  the  values.  In  this 
form  the  information  is  not  readily  useful  for  design.  The  shock  spectrum 
transforms  this  motion  time  listing  into  a  set  of  the  peak  responses  the 
motion  is  able  to  cause  in  a  set  of  simple  mass-spring-dashpot  vibratory 
systems.  Thus,  the  shock  spectrum  describes  the  shock  motion  in  terms 
of  its  capacity  to  excite  simple  vibratory  systems.  This  transformation 
of  the  data  is  now  routinely  accomplished  on  digital  computers;  several 
programs  are  available  in  the  literature  (Ref  3-7). 

However,  the  algebra  required  to  derive  the  calculating  equations 
is  lengthy  and  seldom  published.  For  example,  Lane's  algorithm  (Ref  8), 
which  was  published  in  1964,  is  extremely  efficient,  requires  unbelievably 
few  calculations,  and  is  touiod  as  being  accurate.  It  was  adopted  by 
virtually  all  of  the  aerospace  industry  and  is  still  in  use.  It  was 
derived  by  using  Z  Transform  Theory  (Ref  9)  which  few  understood.  In 
1973,  Cronin  (Ref  4)  published  a  method  of  deriving  Lane's  result  by 
manipulating  a  Duhamel's  integral.  Cronin's  work  showed  the  severity  of 
the  approximation  required  to  obtain  the  Lane  result,  which  raised 
doubts  about  the  accuracy  of  the  method.  O'Hara  (Ref  10),  Nigam  and 
Jennings  (Ref  3),  and  Vernon  (Ref  5)  all  published  two-equation  forms 
for  calculatng  equations.  Their  assumptions  were  clearly  laid  out  and 
satisfactory,  but  computing  with  those  equations  was  comparatively  time 
consuming.  A  comparison  of  Cronin's  work  with  these  two-equation  calcu¬ 
lating  schemes  made  possible  the  development  of  a  new  single-calculating 
equation.  This  new  theory  is  presented  here  for  the  first  time.  The 
new  method  is  more  simply  derived  and  should  contribute  to  a  wider 
understanding  of  the  calculations,  and  thereby  make  them  more  believable 
and  useful. 

COMPUTATIONS 

To  understand  the  computing  to  be  accomplished,  the  stage  must  be 
set.  The  machine  will  be  given  a  list  of  perhaps  1,000  acceleration 
values  in  sequence;  they  au  equally  spaced  values  separated  by  a  known 
time  interval  sampled  from  the  acceleration  versus  time  graph  of  the 
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motion  to  be  analyzed.  The  resulting  spectrum  to  be  calculated  will  be 
a  plot  versus  frequency  of  the  peak  responses  of  simple  vibratory  systems 
of  a  given  damping  ratio  when  exposed  to  that  input  motion.  Figure  1 
contains  two  computed  shock  spectra,  one  for  zero  damping  and  one  for  5% 
damping.  It  is  plotted  on  four-coordinate  paper,  as  discussed  in 
Reference  2.  As  can  be  read  from  the  figure,  this  is  the  shock  spectrum 
of  a  foundation  motion  that  would  cause  an  undamped  10-Hertz  oscillator 
to  attain  a  peak  deflection  of  3  inches,  or  a  100-Hertz  oscillator  to 
attain  a  peak  deflection  of  0.26  inch  and  a  peak  acceleration  of  about 
250g.  Each  curve  on  Figure  1  consists  of  180  equally  spaced  points 
connected  by  straight  lines;  on  that  sized  paper  this  makes  a  very 
smooth  plot.  Now  the  computation  can  be  visualized;  one  must  numerically 
calculate  the  response,  one  at  a  time,  of  180  different  natural  frequency 
vibratory  systems  and  pick  out  and  save  the  peak  value  for  each  frequency. 
A  new  set  of  180  values  is  computed  for  each  damping  ratio.  The  computer 
does  it  quickly  and  inexpensively;  even  the  plotting  is  done  by  machine. 


SINGLE-DEGREE-OF-FREEDOM  SYSTEM  RESPONSE  TO  A  FOUNDATION  MOTION 

The  single-degree-of-freedom  vibratory  system  to  be  conceptually 
used  for  this  data  analysis  is  shown  in  Figure  2.  The  absolute  position 
of  the  foundation  is  y,  and  the  absolute  position  of  the  mass,  m,  is 
given  by  x.  The  spring  stretch,  or  relative  displacement  of  the  mass 
with  respect  to  the  foundation,  is  z,  or 

z  =  x  -  y  (1) 

Dots  are  used  to  indicate  differentiation  with  respect  to  time;  thus,  a 
free  body  diagram  of  the  mass  with  the  reversed  inertia  force  is  as 
shown  in  Figure  3.  Sunning  the  forces  yields 

m  x  +  c(x  -  y)  +  k(x  -  y)  =  0  (2) 
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Writing  Equation  2  in  terns  of  the  relative  displacement,  z,  defined  in 
Equation  1  yields 


mz  +  cz  +  kz  =  -my  (3) 

By  dividing  by  m,  and  using  the  traditional  definitions  and  symbols  in 
Reference  11  given  below 


to  a 

♦  \j  k/m 

(4a) 

C  = 

c 

2  m  to 

(4b) 

t  = 

C/Cc 

(4c) 

one  obtains 

z  ♦  2  £  to  z  +  u)^z  =  -  y  (5) 

Equation  5  gives  the  dynamics  of  the  simple  system  in  terms  of  itB  more 
interpretable  characteristics,  i.e.,  its  undamped  natural  frequency,  u>, 
and  its  damping  ratio,  £.  It  non-sizes  the  equation;  all  single-degree- 
of-freedom  systems  with  the  same  natural  frequency  and  damping  ratio 
must  respond  identically  to  the  same  transient  acceleration.  All  shocks 
or  earthquakes  are  transient  motions,  which  can  be  described  in  terms  of 
their  transient  accelerations.  Thus,  their  effect  on  single-degree-of- 
freedom  systems  is  an  excellent  way  of  organizing  or  classifying  their 
capability  to  damage  equipment. 

The  right  hand  side  of  Equation  5,  the  y,  or  input  to  the  equation, 
will  be  a  list  of  values.  Solutions  for  Equation  5  will  be  developed 
that,  generate  a  list  of  z's  that  result  from  the  input.  The  theory  will 
also  give  us  equations  for  x,  z,  and  z  once  the  list  of  z's  has  been 
found. 

Several  approaches  exist  for  deducing  the  solution  of  Equation  S 
for  a  transient  input,  which  is  what  is  needed  here.  A  good  form,  the 
Duhamel  integral  form,  is  developed  in  most  vibration  texts;  Thomson 
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(Ref  11)  develops  it  for  foundation  excited  motions,  and  O'Hara  (Ref  10) 
gives  it  explicitly  with  initial  conditions  but  without  derivation  as 
follows : 


•  -£u>t 

-Cwt  /  t  .  \  zo  e  . 

2  =  z0  e  p8  wd  1  *  q  8in  wd  t)  *  — ~ —  8in  wd  1 


t 


-  —  J y(x)  sin  u>d(t-T)  dx 


0 

where  z  ,  z  = 

initial  values  of  z  and  z 

o’  0 

md  = 

damped  natural  frequency,  qui 

n  = 

Vi  -  c2 

T  = 

a  dummy  time  variable  of  integration 

(6) 


One  should  not  be  intimidated  by  the  complicated  appearance  of 
Equation  6.  No  one  uses  it  in  this  form  except  to  derive  simpler  rela¬ 
tions.  Equation  6  gives  the  homogeneous  solution  of  Equation  S  in  its 
first  two  terms;  the  motion  described  by  these  terms  is  caused  by  the 
initial  velocity,  zq,  and  displacement,  zqI  existing  at  time  equal  to 
zero.  The  third  term,  the  integral,  is  a  formula  for  finding  the  par¬ 
ticular  solution  for  the  part  of  the  motion)  being  caused  by  the  excita¬ 
tion,  y,  occurring  during  the  time  interval  in  which  Equation  6  is  being 
applied.  The  contribution  to  shock  spectrum  computing  technology  that 
has  been  made  by  this  research  is  the  way  in  which  Equation  6  was  used. 
In  the  next  few  paragraphs,  two  very  popular  computing  equations  are 
reviewed.  All  of  the  results  presented  below  are  derived  and  further 
discussed  in  Appendix  A. 


DISPLACEMENT  SINGLE  RECURSIVE  EQUATION,  TRAPEZOIDAL  RULE  INTEGRAL 
APPROXIMATION 

The  simplest  of  the  calculating  solutions  will  be  discussed  first, 
and  will  serve  to  illustrate  the  use  of  all  the  calculating  solutions. 
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It  was  first  presented  by  Lane  (Ref  8)  and  is  in  current  use  by  most  of 
the  aerospace  community  (Ref  4,12,13).  The  result  is  given  below 


z . 


X 


z--2 +  ai6  Vi +  *i  y 


15 


i-1 


(7) 


where  the  constants  are  given  by 


15 


"16 


2  e  cos  u).  h 
a 


*1  = 


.  -£iuh  . 

he*  sin  uij  h 

w. 


(7a) 

(7b) 

(7c) 


and  where  h  is  the  time  interval  between  samples  of  y.  This  is  derived 
in  Appendix  A  as  Equation  A- 11a. 

The  use  of  Equation  7  can  be  envisioned  as  follows.  Consider  that 
for  a  given  damping  value,  one  wants  to  compute  the  peak  z  for  some 
frequency  u»j.  The  constants  a^,  a^,  and  Yj  are  first  calculated. 

Then  y^  is  substituted  into  Equation  7  with  z^  and  zq  equaling  zero. 

This  yields  a  value  for  Next  y^  and  z ^  are  used  to  compute  z^  with 

Zj  again  taken  equal  to  zero.  The  value  z^  is  calculated  from  Equation  7 
by  using  y^,  z^»  and  z^;  and,  thus,  the  process  of  calculating  the  z's 
from  the  y's  continues.  At  each  step  of  the  computation,  each  newly 
computed  value  of  z  is  tested  to  see  if  it  should  replace  what  up  until 
that  time  has  been  the  most  positive  and  negative  values.  The  absolute 
largest  of  the  two  is  the  maximum.  When  the  processing  of  the  list  of 
y's  has  been  completed  and  the  maximum  z  found,  one  has  found  the  maximum 
"during"  values  of  z  for  that  frequency.  This  process  is  repeated  for 
each  frequency  at  which  a  shock  spectrum  value  is  desired. 

One  can  refer  to  Equation  7  as  a  displacement  recursive  equation 
for  use  with  equally  spaced  digitized  values  of  acceleration.  The  list 
of  displacements  computed  therefrom  completely  defines  the  resulting 
motion.  In  Appendix  A  it  is  shown  that  if  the  velocity  at  any  time  is 
desired,  it  can  be  computed  from  these  displacement  values  as  follows 
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®11  2i-l  *  ®12 


♦  Y 


2  yi 


(7d) 


where  Che  constants  are  given  by 

au  =  -(n  w  e"^“*l)/(sin  qwh) 

a  *  u»(q/tan  rnuto  *  C) 

12 

Y2  =  -  h/2 

If  the  relative  acceleration,  z^,  is  of  interest,  it  is  obtained  by 
using  the  values  z^  and  y^,  along  with  the  computed  value  of  z^  in 
Equation  5,  which  gives 

zi  =  -  2  C  u.  -  ui2  z£  -  y .  (7e) 


If  the  absolute  acceleration  is  desired,  one  notes  from  Equation  1  that 


=  zi  *  yi 


(7f) 


and,  thus,  from  Equation  7e  one  finds 


=  •  2  t  w 


u> 


(7g) 


Equations  7,  the  aerospace  industry  equations,  are  a  complete  set. 

References  4  and  13  give  Fortran  program  lists  for  their  use.  As  men* 
tioned  in  Appendix  A,  the  apparent  crudeness  of  the  approximate  Integra* 
tion  method  was  not  made  clear  until  Reference  4  was  published  in  1973. 
Nothing  has  been  written  criticizing  these  equations,  including  Reference  4, 
except  that  they  are  completely  ignored  by  the  respected  earthquske 
analysis  community  (Ref  3  and  7).  Reed  (Ref  13)  compares  several  com* 
puting  approaches  and  finds  these  equations  inaccurate  for  coarse  sampling 
rates,  but  otherwise  acceptable.  The  computational  approach  advocated 
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here  resulted  from  trying  to  reconcile  this  computing  approach  with  that 
of  References  3,  7,  and  10.  A  suitably  simple  theory  was  found  for 
transforming  the  conventional  equations  (Ref  3  and  10)  into  a  single 
displacement  recursive  equation  so  that  simplicity  of  a  single  equation 
could  be  retained  along  with  the  accuracy  of  the  more  conventional 
integration  approximations . 


DISPLACEMENT  AND  VELOCITY  PAIR  OF  RECURSIVE  EQUATIONS,  ACCELERATION 
APPROXIMATED  AS  A  STRAIGHT  LINE  IN  INTEGRAL 


The  more  conventional  equations  referred  to  above  were  first  pre¬ 
sented  by  O'Hara  (Ref  10)  in  1962  and  then  independently  again  by  Nigaa 
and  Jennings  (Ref  3)  in  1968.  They  approximate  the  acceleration  by  a 
straight  line  between  the  sample  points  ^nd  then  integrate  the  Integral 
of  Equation  6.  This  yields  a  pair  of  equations  for  processing  the  data. 
They  are  derived  in  Appendix  A  as  Equations  A-18c  and  A-I8d  and  are  as 
follows 


Zi  *  “l  2i-l  +  °2  h-l  +  °29  *1-1  +  “30  h 
\  =  °4  Zi-1  +  “5  *1-1  +  °31  h-l  *  a32  *1 


(8a) 

(8b) 


The  constants,  which  are  a  complicated  expression,  but  of  the  same 
variables  as  in  Equation  7,  are  given  in  Appendix  A  as  Equations  A-4c, 
A-4d,  A-4f,  A-4g,  A-181,  A-18J,  A-I8k,  A-181.  In  this  computing  method, 
two  outputs,  both  z  and  z,  must  be  computed  for  each  step,  and  then  used 
in  the  next.  Since  both  z  and  z  are  being  constantly  computed,  it  is  a 
simple  matter  to  use  Equation  7g  to  compute  the  acceleration  of  the 
mass,  ft,  at  each  point  if  desired.  This  is  the  procedure  currently  in 
use  for  most  of  the  U.S.  earthquake  data  (Ref  7),  and  hence  is  a  highly 
respected  spproacb.  Note  that  for  each  step,  or  for  each  point  of  input 
data,  eight  multiplications  and  six  additions  are  required.  What  was 
found  during  the  course  of  this  work  was  that  the  above  pair  could  be 
exactly  reduced  to  a  single  equation;  this  results  in  each  step  only 
requiring  five  multiplications  and  four  additions. 
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DISPLACEMENT  SINGLE  RECURSIVE  EQUATION,  ACCELERATION  APPROXIMATED  AS  A 
STRAIGHT  LINE  IN  INTEGRAL 


la  Appendix  A  it  is  shown  thst  with  some  shifting  of  the  indices, 
one  csn  combine  Equations  8a  and  8b  into  the  single  equation  given  below 


*15  Zi-2  +  “16  *i-l  +  a24  *i-2  +  °25  Vl  +  a26  *i 


(9a) 


The  constants  again  are  complicated,  but  easily  computed  and  are  given 
in  Appendix  A  as  Equations  A-lOe,  A-lOf,  A-19b,  A-19c,  and  A-19d.  Note 
that  the  constants  applied  to  the  z's  are  the  same  as  those  of  Equation  7, 
Lane's  result  (Ref  8).  No  matter  how  sophisticated  an  approximate 
integration  is,  these  same  two  constants  appear.  It  does  seem  more 
reasonable  that  three  input  values  are  required  for  each  step,  including 
the  value,  y^,  for  the  time  at  which  the  output,  z^,  is  being  computed. 

As  was  the  case  previously,  the  data  or  list  of  y^'s  is  marched 
through  Equation  9a,  which  yields  a  z±  for  each  y^  If  the  velocity  at 
any  point  is  desired,  it  can  be  computed  as  shown  in  Appendix  A, 

Equation  A~20a,  from  the  z's  as  follows 

*1  *  °11  *1-1  ♦  *12  *1  *  *27  n-1  *  *28  h  <,b) 

The  constants  are  given  in  Appendix  A  as  Equations  A-lOa,  A-lOc,  A-20b, 
and  A-20c.  If  the  absolute  acceleration  of  the  mass,  x^,  is  desired,  z^ 
is  computed  and  used  with  z^  in  Equation  7g. 

RESIDUAL  SPECTRUM  CALCULATION 

As  has  been  discussed,  the  spectrum  computation  consists  of  exciting 
the  theoretical  model  of  Figure  1,  which  is  described  by  the  differential 
Equation  5,  by  means  of  a  calculation  algorithm  (such  as  Equations  7,  6, 
or  9)  with  the  excitation  given  as  a  sequence  of  equally  spaced  values. 
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After  the  excitation  is  over,  the  oscillator  continues  to  ring,  and  the 
low  frequency  oscillators  attain  their  greatest  values  after  the  excita¬ 
tion  has  passed.  In  the  previously  discussed  algorithms,  the  calculation 
of  these  "residual  values"  is  accomplished  by  continuing  the  algorithm 
with  zero  excitation  values  for  one  period  of  the  oscillator.  However, 
this  is  not  the  only  way  of  doing  it,  or  necessarily  the  best  way. 

Often  one  computes  spectral  values  for  oscillators  with  periods  far 
longer  than  the  excitation  duration,  and,  in  this  case,  the  residual 
integration  region  is  necessarily  longer  than  the  original  pulse. 

• 

Another  approach  is  to  compute  the  final  velocity  and  displacement,  z 
and  z,  at  the  end  of  the  excitation,  and  then  calculate  the  positive 
maximum  and  negative  minimum  from  the  free  vibration  solution.  This  is 
the  procedure  used  in  the  new  shock  spectrum  computing  method  and,  thus, 
will  be  derived  here  with  some  of  the  detail  given  in  Appendix  C. 

The  solution  for  a  damped  oscillator  undergoing  free  (or  decaying) 
vibrations  can  be  taken  from  the  first  two  terms  of  Equation  6  (since  y 
is  then  zero) ,  or 

z  =  z  e  ^jcos  10 .  t  +  4-  sin  u>.  t) 
o  \  d  q  d  f 

-£« ut 

2  e 

+  ———sin  u).  t  (10) 

u».  d 


By  rearranging,  this  can  be  put  in  a  form  more  convenient  for  manipulating 
as  follows: 


z 


e-C  wt 


(A  sin  u> .  t  ♦  B  cos  u>.  t) 
a  a 


where  A 


B  = 


z 

0 


(11) 


(11a) 

(lib) 


It  can  be  shown  by  substitution  that  a  convenient  form  for  the  derivative 

a 

of  Equation  11  or  z  is 
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z 


(12) 


=  u>  e  ^Ult  (A  cotCiOj  t  +  6)  -  B  sin(u^  t  +  6)) 


where  sin  6  = 

(12a) 

cos  6  = 

n 

(12b) 

0  = 

Vi  -  t2 

(12c) 

u>d  = 

n  ui 

(12d) 

Since  Equation  11  is  a  damped  vibration  or  at  the  most  with  £  =  0,  a 
constant  vibration,  both  the  nost  positive  and  the  most  negative  value 
must  occur  in  the  first  cycle  or  at  t  =  0.  Thus,  z  and  z  .  must  be 

pursued  in  the  first  cycle.  If  not,  when  t  *  0,  these  occur  when  z 
equals  zero  or  when 

A  coa(u>d  t  +  6)  =  B  sinCui^  t  +  6)  (13) 

If 

0  =  u»d  t  +  6  (14) 

one  seeks  0  such  that 

A  cos  0  *  B  sin  0  (15) 

or 

tan  0  *  A/B  (15a) 

and  also  seeks  two  values  0j  and  0^  such  that  the  results  (u>dt)j  and 
(u»dt)2  are  between  zero  and  2 n.  Clearly,  0j  and  02  are  consecutive 
angles  at  which  the  velocity  goes  to  zero;  thus  (since  the  tangent  goes 
through  zero  every  n  radians), 


?!  +  n 


(16) 


P 


2 


Now  one  oust  go  through  an  amount  of  detail  to  assure  the  finding  of  the 
very  first  time  after  time  equals  zero  that  z  equals  zero.  Nine  cases 
must  be  considered  as  indicated: 


This  is  done  carefully  in  Appendix  C.  The  results  are  given  below  as 
the  set  of  "if"  statements  used  in  programming  the  computation. 

1.  If  both  A  and  B  are  zero,  no  residual  response  results. 


2.  If  not,  and  A  *  0,  then  Pj  =  0.  (17a) 

3.  If  not,  and  B  =  0,  then  pj  -  ■j.  (17b) 

A.  If  not,  and  A  and  B  have  the  same  sign, 

=  Tan-1  (A/B)  (17c) 

S.  If  not, 

pj  =  n  -  Tan"1  (-A/B)  (17d) 


Now  with  this  value  of  P^,  one  final  test  is  required.  If,  and  only  if 
P  -  <?  <  0  (17e) 

one  must  add  n  to  P^,  or 

p,  -  <e,>ou  *  *  <«« 
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Now  with  the  acceptable  value  of  Pj,  (w^t) ,  can  be  obtained  with  Equation  14 
as  follows 

(WjOj  5  Pj  -  5  U7g) 

The  second  value  of  (ii^t)  is  obtained  from 

(u»dt)2  =  (WjOj  +  n  (17b) 

These  values  of  (w^t)^  and  (WjOj  are  substituted  successively  in*.o 

Equation  11,  which  will  yield  values  for  z  and  z  .  .  These  values 

max  Bin 

are  compared  with  z  to  make  sure  it  is  not  greater  than  2  or  less 

o  max 

than  za±at  And,  thus,  the  residual  spectrum  values  are  computed.  Sub¬ 
routine  RESID  listed  at  the  end  of  Program  SPCTRM  in  Appendix  B  uses 
this  procedure  to  compute  the  residual  spectrum  values. 

Some  additional  explanation  of  the  need  for  the  residual  spectrum 
computation  procedure  is  given  here,  since  the  previously  mentioned 
programs  do  not  go  to  this  trouble.  To  get  the  residual  response,  one 
must  first  use  a  numerical  method  (such  as  Equations  7,  8,  or  9)  to  find 
the  displacement  and  velocity  (zq  and  zq)  at  the  end  of  the  "during" 
portion  of  the  response.  The  other  computations  continue  the  numerical 
procedure  to  find  the  residual  maximum  values.  This  procedure  uses  a 
theoretically  exact  method  to  compute  the  values;  it  is  faster,  more 
accurate,  but  a  little  clumsier  to  program. 


THE  COMPUTER  PROGRAMS 

Appendix  B  gives  the  Fortran  listings  for  the  two  computer  programs 
used  in  the  preparation  of  shock  and  response  spectra  by  this  new  method. 
The  first,  SPCTRM,  computes  the  shock  spectrum  and  lists  as  output  the 
values  of  the  maximum  and  minimum  displacement  and  pseudovelocity  for 
each  frequency  considered.  The  second,  PLTVLF,  is  the  program  used  to 
plot  the  maximum,  positive  or  negative  shock  spectrum  on  the  four- 
coordinate  paper.  Note  that  the  precise  frequencies  used  for  calculation 
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*  >»WifT 


are  selected  by  the  program  so  that  they  will  be  equally  spaced  when 
plotted  logarithmically.  Appendix  D  documents  the  logic  for  this  section 
of  program  SPCTRM.  Otherwise  the  steps  and  symbols  of  the  program 
coincide  with  the  report  body  and  appendices,  and  should  be  completely 
readable.  The  programs  are  in  current  use  on  a  CDC  time-sharing  terminal; 
the  plotting  is  done  on  an  interconnected  Hewlett  Packard  Model  7202A 
Graphic  Platten.  The  format  statements  show  the  form  of  the  input  data. 
The  programs  interrogate  the  user  and,  thus,  supply  the  needed  documen¬ 
tation. 


COMMENTS 


The  Accuracy  of  the  Algorithms 


Equation  9a  is  more  accurate  than  any  others  thus  far  published. 
Because  these  recursive  calculating  equations  use  previously  computed 
results  for  each  succeeding  step,  they  propagate  any  errors  introduced. 
Thus,  if  fewer  multiplications  and  additions  are  required  for  each  step, 
there  is  less  opportunity  for  error  accumulation  due  to  roundoff  and 
subtraction  of  nearly  equal  values. 

Of  the  algorithms  reviewed,  Equation  7,  the  "trapezoidal"  approach, 
has  the  fewest  multiplications  and  additions.  However,  it  assumes  each 
factor  of  the  integral  is  constant  during  the  integration.  If  one 
graphs  the  value  of  the  factors  over  the  integration  interval  and  then 
compares  an  estimate  of  the  integral  of  their  product  with  that  assumed 
by  the  algorithm,  the  crudeness  of  the  approximation  becomes  apparent. 
Thus,  the  "trapezoidal"  rule  approach  was  discarded  because  of  these 
initial  assumptions . 

Conversely,  if  one  graphs  the  values  of  the  factors  similarly  with 
the  straight  line  integration  approach,  one  can  infer  a  convergence  to 
the  "true"  value  as  the  steps  become  smaller.  Such  an  examination  even 
leads  one  to  assume  the  straight  line  approach  is  exact;  such  a  comment 
was  even  made  in  Reference  3.  That  is  why  most  of  the  algorithms  use 
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the  straight  liae  approach.  Of  these  the  new  single  recursive  equation 
developed  here  has  by  far  the  fewest  multiplications  and  additions  and, 
thus,  is  the  most  accurate. 

A  rigorous  study  of  the  accuracy  of  the  approaches  has  not  been 
undertaken,  and  would  certainly  be  welcomed.  However,  perfunctory 
accuracy  statements  that  pretend  to  bound  the  error  by  referring  to 
higher  derivatives  of  the  input  are  not  applicable  to  digitized  data. 

Nor  is  there  any  way  in  which  one  can  construe  the  data  to  be  composed 
of  low  order  polynomials.  As  is  discussed  further  in  Appendix  A,  one 
does  not  want  to  sample  any  coarser  than  about  10  points  per  cycle  of 
the  highest  frequency  present.  At  this  high  a  sampling  rate,  the  straight* 
1 ine-be tween- the -point b  assumption  makes  an  accurate  appearing  plot  of 
the  data,  which  leaves  one  to  believe  the  integration  has  been  adequately 
approximated. 

Allowable  Frequency  Ranges 

This  new  program,  as  well  as  all  other  shock  spectrum  computer 
programs,  has  limitations  on  the  frequencies  for  which  it  can  compute 
accurate  shock  spectrum  values.  The  allowable  frequencies  are  dependent 
upon  the  sampling  rate  of  the  data  used  as  input  and  upon  the  number  of 
digits  the  particular  computer  uses  in  performing  its  arithmetic  opera¬ 
tions.  The  high  frequency  limit  for  the  programs  discussed  in  this 
report  is  one  tenth  of  the  sampling  rate.  This  assures  a  5%  accuracy 
and  is  discussed  by  most  authors  (Ref  3). 

That  a  low  frequency  limit  exists  has  not  been  published  up  to  now. 

The  constants  that  apply  to  the  input  data  in  Equations  7,  8,  and  9a 
contain  sines  and  cosines  with  arguments  that  become  very  small  at  low 
frequencies  and  lead  to  the  subtraction  of  nearly  equal  numbers  (e.g., 
equal  to  the  12th  or  13th  digit  in  a  14-digit  number).  At  very  low 
frequencies  this  can  lead  to  the  computer  generating  ridiculous  values 
off  by  a  factor  of  ten  or  more.  Time  did  not.  permit  a  thorough  investi¬ 
gation  of  this  effect,  but  an  empirical  testing  of  the  constants  did 
indicate  that  for  the  computer  used,  which  calculates  to  14-digit  accuracy 
in  single  precision,  accuracy  could  be  maintained  for  frequencies  as  low 
as  one  eight- thousandth  of  the  sampling  rate. 
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These  frequency  limitations  are  more  completely  discussed  in 
Appendix  A.  In  symbols,  the  range  of  allowable  frequencies  for  the 
programs  can  be  expressed  as 

f 

10  S  ~  S  8000  (18) 

Thus,  for  data  that  have  been  sampled  at  10,000  points  per  second,  one 
could  compute  shock  spectrum  values  for  frequencies  as  high  as  1,000  Hertz 
and  as  low  as  1.25  Hertz.  The  current  program  does  not  impose  this 
restriction;  it  is  the  responsibility  of  the  user. 


SUMMARY 

The  body  of  the  report  has  presented  a  background  in  shock  spectrum 
calculation  methods  applicable  to  digital  computation  that  has  not  been 
available.  In  general  terms  in  the  text  and  in  detail  in  Appendix  A, 
the  two  most  popular  methods  and  a  new  method  are  described,  defined, 
and  completely  derived.  This  new  method  decreases  the  numbers  of  addi¬ 
tions  and  multiplications  required  by  40%  and,  thereby,  increases  accuracy 
while  decreasing  cost.  The  main  purpose  of  the  report  is  to  document 
this  new  method  and  make  it  available.  Appendix  B  gives  computer  program 
lists  for  this  new  method. 

There  are  still  unresolved  questions  and  certainly  room  for  improve¬ 
ment;  these  matters  are  discussed.  In  particular,  a  new  potential 
source  of  inaccuracy  was  discovered  in  trying  to  extend  to  earthquake 
frequencies  data  that  had  been  sampled  for  explosion  frequencies.  The 
only  way  one  can  compute  spectra  for  the  very  low  earthquake  frequencies 
is  to  have  the  data  digitized  at  an  appropriate  sampling  rate. 

The  author  believes  that  a  greater  understanding  of  the  shock 
spectrum  is  valuable  for  evaluating  equipment  installation  in  dynamic 
environments.  Mere  plotting  of  the  shock  spectrum  on  four-coordinate 
paper  gives  a  vastly  improved  appreciation  of  the  damaging  capacity  of  a 
given  environment.  Required  rattlespace  stands  out  and,  to  the  extent 
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that  one  can  approximate  the  gross  natural  frequency  of  the  equipment  on 
its  support,  the  peak  breakaway  forces  can  be  estimated.  A  considerable 
increase  in  safety  and  facility  hardness  from  the  effects  of  explosion 
and  earthquake  is  available  inexpensively  from  increased  understanding 
and  use  of  the  shock  spectrum. 

Finally,  this  study  of  the  shock  spectrum  analysis  of  sampled 
acceleration  time  histories  not  only  reduced  and  clarified  the  required 
calculations,  it  also  provoked  several  ideas  for  even  further  reducing 
the  program  size.  It  must  be  assumed  that  soon  small  calculators  will 
be  doing  the  analysis.  Inexpensive  digitizers  and  plotters  will  follow, 
and  all  will  have  access  to  many  new  toola  for  intelligent  use  of  the 
data  that  now  is  only  plotted. 
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Appendix  A 


DERIVATION  OF  RECURSIVE  EQUATIONS 


DERIVATION  APPROACH  AND  NEW  SYMBOLS 

As  he*  been  discussed  in  the  text,  one  seeks  s  solution  of  Equstion  5, 

z  ♦  2  (  ut  z  +  ui^  z  *  -y  (5) 

spplicsble  to  s  list  of  closely  snd  equslly  spsced  vslues  of  the  input 
foundstion  sccelerstion,  y.  Hence,  three  such  solutions  of  Equstion  5 
will  be  developed  that  apply  from  one  excitation  point  to  the  next. 

They  will  be  recursive  calculating  equations,  through  which  the  input 
dsts  sre  successively  psssed,  thereby  yielding  s  list  of  output  values, 
z.  The  text  equations  snd  definitions  1  through  6  will  be  used.  The 
approach  of  this  Appendix  is  to  give  s  general  beginning  leading  to  two 
known  results  snd  the  new  result. 

The  Duhamel  integral  solution,  Equation  6,  is  the  point  of  departure. 
The  derivative  of  this  solution  with  respect  to  tine  is  needed,  and  this 
oust  be  done  carefully  since  t  is  a  Unit  of,  and  occurs  within,  the 
integral.  Pipes  (Ref  14)  explains  this  carefully,  snd  following  bis 
instruction  one  obtains 


z  *  -  —2 —  e  sin  ut .  t  ♦  ~  e’^^-C  sin  u>.  t 

n  d  n  *  d 


v 

+  f|  cos  u>d  t)  -  — ■ yV(t)  e"^w^t"x^J-C  sin  w^t-t) 


+  n  cos  WjCt-t))  di 


(A-l) 
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naomm  mm  a*m-mr  mm d 


Consider  applying  Equations  6  and  A-l  from  one  value  of  y  to  the  next; 
i.e.,  from  y^  to  y^.  As  can  be  read,  the  equations  presume  zq  and  zq 
are  known.  Equation  6  would  give  z^,  and  Equation  A-l  would  give  z^. 

To  clarify  what  follows,  some  new  symbols  must  be  defined.  The 
time  interval  between  the  equally  spaced  values  will  be  h.  Equations  6 
and  A-l  are  to  be  applied  when  t  equals  h.  Therefore,  define 


it*  s  e  sin  u>d  h 

»  e  cos  u»d  h 

h 

01  s  J ir(x)  e"t“*Ol"X)  gla  ^(h-x)  dt 

o 

h 

01  b  y* y(x)  e  ^h'T^  cos  uid(h-t)  dx 
0 


(A-2a) 

(A-2b) 


(A-2c) 


(A-2d) 


Using  these  new  symbols,  Equations  6  and  A-l  become 

(A- 3a) 
(A- 3b) 

Note  that  as  these  equations  sre  applied  to  the  input  list  of  y's, 
h,  w,  Ci  u>d,  OfX,  and  i)»  are  constants.  If  subscripted  a's  are  to 
represent  constants,  Equations  A-3  become 


*1  *  “l  ♦  °2  *  “3  S01 


z.  =  a.  z  +  a,  z  +  a,  SA,  +  a,  Cn, 
1  4o  5o  6  01  7  01 


where  dj  s  x  +  ^  <|< 


(A-4s) 

(A-4b) 

(A-4c) 
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jl. 

(A-4d) 

n  uj 

-  i 

(A»4e) 

n  w 

U)  t|l 

"  ~n~ 

(A-4f ) 

-11  +  x 
n 

(A-4g) 

t 

n 

(A-4h) 

-l 

(A-4i) 

Consider  the  integrals  5  and  C-  The  way  in  which  they  are  evaluated 
constitutes  the  differences  between  the  various  results. 


TRAPEZOIDAL  APPROXIMATION 

The  crudest  method  in  popular  use  is  the  trapezoidal  rule  (Ref  4), 
which  states  that  the  integrand  is  taken  to  be  constant  over  the  interval 
and  equal  to  the  average  of  the  values  of  the  integrand  at  the  beginning 
and  end  of  the  interval. 

From  Equation  A-2c, 

h 

S01  =  f  ^(t)  e"kU(h_T)  si"  «>d(b-T)  (A-2c) 

0 

The  average  of  the  end  point  values  is 

Aj  =  -j-  ^yQ  e”^ul*1  sin  u>dh  ♦  0^ 
or 

Aa  =  -|y0  (A-5a) 
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Using  this  as  the  integrand  in  A-2c  yields 


(A-5b) 


for  the  value  of  Sqj  by  the  trapezoidal  rule.  To  get  the  value  of 
Cqj,  consider 

h 

CQ1  =  cos  tud(h-T)  dx  (A-2d) 

0 


The  average  of  the  end  point  values  is 

A2  =  i  (?o  e"twh  C08  wd  h  +  *l) 

Using  this  as  the  integrand  in  A-2d  yields 
coi  ■ 


(A-5c) 


(A-5d) 


The  point  to  note  here  is  that  in  this  case,  and  in  the  other  cases 
to  be  considered,  the  integrals,  SQ1  and  CQ1,  are  composed  of  constants 
and  input  data  values;  no  z's  or  output  values.  This  can  be  seen  more 
clearly  as  follows. 

By  substituting  these  values  of  and  Cqj,  Equations  A-5b  and 
A-5d  into  Equations  A-4a  and  A-4b,  one  obtains 


Z1  =  al  zo  *  °2  zo  +  °3  *0  (A- 6a) 

and 

Z1  =  °4  *o  +  a5  zo  +  a6  ^0  +  °7  ^1  (A-6b) 

where  these  new  a's  are  defined  to  be 
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(A-6c) 


Oj  S  <*3  «|(  h/2 

°i  2  T  («6  *  ♦  “7) 

a}  s  or7  h/2 


(A-6d) 
(A-6 e) 


Thus,  it  is  seen  that,  upon  evaluation  of  Sg7  and  Cq7>  which  become 
products  of  constants  and  the  input  y's,  Equations  A-4a  and  A-4b  become 
Equations  A-6a  and  A-6b.  These  form  a  pair  that  are  used  together  to 
generate  a  list  of  both  z's  and  z's,  from  the  y's.  Assuming  the  fit's 
have  been  computed  and  the  full  list  of  the  y's  is  on  hand,  one  would 
start  the  calculating  at  yQ,  taking  zq  and  zq  equal  to  zero,  and  proceed 
to  z^.  Knowing  and  z^  the  pair  of  equations  will  yield  zi+1  and 
zi+1  and  so  on. 

With  this  in  mind  it  is  now  possible  to  go  back  and  combine  the  two 
equations  (A-4a  and  A-4b)  into  a  single  equation  for  z.  alone,  thereby 
eliminating  the  need  to  keep  computing  both  z  and  z.  This  development 
does  not  require  Z  Transform  Theory  (Ref  8  and  9)  and  is  more  straight¬ 
forward  than  the  approach  of  Cronin  (Ref  4} . 


DEVELOPMENT  OF  A  SINGLE  RECURSIVE  EQUATION 

One  begins  by  solving  Equation  A-4a  for  z  and  defining  new  cr's. 
This  yields 


*.  *  "»  zo  *  “9  21  *  °10  S01 


where  a 


8 


(A-7a) 


(A-7b) 

(A-7c) 


(A-7d) 
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Equation  A-7a  is  next  substituted  into  A-4b,  which  then  becomes 


* 1  1  °4  *.  *  “5(<,8  *o  *  “s  *1  *  “10  S01> 

*  “6  S01  *  “?  C01  <*-*•> 

By  defining  new  a's,  this  can  be  written  as 


*1  s  all 

z 

0 

+  aJ2  Zj  +  al3  SQ1  ♦  aJ4  CQJ 

(A-8b) 

where  = 

*4 

*  a5  *8 

(A-8c) 

a12  2 

°5 

a9 

(A-8d) 

a13  2 

a5 

“lO  +  a6 

(A-8e) 

°14  2 

°7 

(A-8f) 

This  is  a  very  interesting  result.  It  states  that  the  velocity  at  point 
1  is  exactly  given  by  zQ  and  z^,  and  the  values  of  the  two  integrals. 
Note  that  it  does  not  require  zq. 

The  next  step  is  to  realize  that  Equation  A-4a  also  applies  between 
points  1  and  2  equally  as  well  as  between  points  0  and  1,  and,  thus,  the 
following  can  be  written 


•l  *1  *  °2  *1  *  "3  S12 


(A-9a) 


The  value  of  Zj  from  Equation  A-8b  is  substituted  into  this  form,  which 
yields 
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And  finally  by  defining  sore  a's,  this  becomes 


Z2  =  “lS 

zo 

+  al6  Z1  *  *17  S01  +  “18  S12  *  al9  C01 

(A-9c) 

where  aJ5  = 

a2 

“il 

(A-9d) 

Q 

O' 

III 

al 

+  <*2  aJ2 

(A-9e) 

Q 

Nl 

III 

a2 

a13 

(A-9f) 

a!8  5 

a3 

(A-9g) 

°19  2 

a2 

al4 

(A-9h) 

Equation  A-9c  is  the  new  result.  Note  that  it  is  quite  different 
from  Equation  A-4a  in  that  it  requires  no  input  velocity;  thus, 

Equation  A-4b  does  not  have  to  be  computed  along  with  it.  It  reduces 
the  number  of  computations  and,  thus,  the  accumulation  of  error.  It  is 
therefore  more  accurate  and  economical.  If  the  velocity  is  desired,  it 
can  be  calculated  with  Equation  A-8b.  As  it  stands  it  is  exact.  Approx¬ 
imation  will  be  introduced  when  the  integrals  are  evaluated  and  by 


round-off  error  accumulated  in  the  calculations. 

Before  using  Equation  A-9c  with  various  approximations  for 
values  of  the  integrals,  and  Cqj,  many  of  the  or  values  will 
required  and  so  the  most  important  ones  are  listed  below. 

the 

be 

all 

= 

f|  u»  _-2frjuh  _  n«»e 

«j»  sin  u»j  h 

(A- 10a) 

“12 

= 

“  t  *  “  Wu-E  •  t 

(A- 10b) 

a13 

r 

X  _  1 

4*  tan  Wj  h 

(A-lOc) 

al4 

= 

-1 

(A-lOd) 

al5 

= 

CM 

1 

t> 

1 

(A-lOe) 
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“l6 

=  2  X 

(A-10f) 

“17 

_  X 
n  ut 

(A-lOg) 

“13 

_  -  i 

n  u> 

(A-lOh) 

“19 

n 

(A-lOi) 

TRAPEZOIDAL  RULE  RECURSIVE  EQUATIONS 

The  calculating  equations  that  result  from  approximating  the  integrals 
by  the  trapezoidal  rule  are  obtained  by  substituting  Equations  A-5b  and 
A-Sd  into  Equation  A-9c.  This  yields 

*2  *  “15  z»  *  “16  21  *  h  »1 

where 

Yx  =  -(h  (|<)/(0  w)  (A-llb) 

This  is  the  result  given  in  References  4  and  8.  Lane  (Ref  8)  derived  it 

in  1964  by  Z  Transform  Theory,  which  did  not  make  clear  the  approximation 
used  in  evaluating  the  integral.  Cronin  (Ref  4)  in  1973  was  the  first 
to  derive  this  result  by  conventional  means;  his  work  showed  that  it 
could  be  derived  with  the  trapezoidal  rule.  It  is  indeed  simple:  from 
only  one  input  value,  and  not  even  the  value  of  the  input  corresponding 
to  the  time  at  which  the  output  is  being  computed,  and  two  previous 
output  values,  the  equation  yields  a  new  output  value. 

The  velocity  equation  is  given  by  substituting  Equations  A-5b  and 
A-5d  into  Equation  A-8b.  This  yields 

21  =  “ll  2o  +  “l2  Z1  +  h  *1  (A-llc) 

where 

Y2  =  -  h/2  (A-lld) 
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This  equation  is  used  to  compute  a  velocity  if  requited;  e.g.,  the  final 
velocity  is  required  as  input  to  the  residual  response  calculation. 
Reiterating,  Equation  A-lla  is  used  to  compute  a  list  of  z's  from  the 
list  of  y's;  i.e.,  if  y^  is  the  first  non-zero  point,  Equation  A-lla 
yields  z^  as  the  first  non-zero  output  value.  Next  y^  and  z^  are  uaed 
to  compute  z^;  next  y^,  z^,  and  z^  are  used  to  compute  z^,  etc. 

Equation  A-lla  written  for  the  general  (n+1)  point  is 

Vl  *  ”15  Vl  *  ”16  *«  *  h  \  (A-U”> 

which  is  a  more  common  way  of  writing  a  recursive  equation. 


STRAIGHT  LINE  ACCELERATION  APPROXIMATION  IN  INTEGRAL 

There  are  several  other  approaches  one  could  take  in  evaluating  the 
integrals  S  and  C.  Cronin  (Ref  4)  suggested  considering  the  values  of 
the  integrand  at  three  consecutive  points  as  defining  a  parabola  and 
then  integrating  the  parabola  so  defined.  Rather  than  approximate  the 
whole  integrand  by  an  approximating  function,  O'Hara  (Ref  10)  in  1962, 
and  Nigam  and  Jennings  (Ref  3)  in  1968,  both  derived  the  pair  of  computing 
equations  that  result  from  Equations  A-4a  and  A-4b  when  the  acceleration 
input  is  approximated  as  a  straight  line  between  two  adjacent  values. 

They  then  integrated  the  product  of  the  approximate  acceleration  function 
and  the  remainder  of  the  integrand.  Vernon  (Ref  5)  in  1967  published  an 
unusual  pair,  similar  to  Equations  A-4a  and  A-4b,  again  with  the  straight 

line  acceleration  approximation,  except  he  chose  to  compute  x  and  x  as 

• 

output  values  (as  opposed  to  z  and  z).  Rodeman  (Ref  15)  in  1974  used  Z 
Transform  Theory  more  accurately  than  Lane  (Ref  8),  approximating  only 
the  acceleration  as  a  straight  line,  and  deriving  a  single  recursive 
equation  with  acceleration,  x,  as  the  output,  aB  opposed  to  z.  Thus, 
the  straight  line  approximation  has  a  host  of  advocates.  Except  for 
Cronin  (Ref  4),  none  found  the  single  recursive  equation  without  Z 
Transform  Theory. 
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A*  to  the  accuracy  of  the  straight  line  approximation,  Reference  3 
calls  the  method  "exact".  O'Hara  (Ref  10)  gives  formulas  for  approxi¬ 
mating  the  acceleration  by  a  parabola,  which  clearly  has  to  be  better 
than  a  straight  line.  The  philosophy  taken  here  is  that,  if  the  sequence 
of  digits  when  plotted  with  straight  lines  connecting  adjacent  points  is 
a  good  visual  approximation  to  the  analog  data,  then  the  digitizing  is 
adequate  for  approxiatating  the  acceleration  function  with  straight 
lines . 

To  evaluate  the  integrals  Sqj  and  Cqj  with  the  acceleration  taken 
as  a  straight  line  between  adjacent  values,  or  yQ  and  y  ,  let 

y(t)  =  a  ♦  b  x  (A-12a) 

where 

*  s  y0  (A-12b) 

b  s  (?!  -  y'0)/h  (A-12c) 

Equation  A-12a  is  used  in  the  definition  of  S01>  Equation  A-2c,  to  yield 
h 

=  J" (a  +  b  x)  ain  WjCh-x)  dx  (A-13a) 

0 

To  integrate  this  expression,  change  the  variable 
u  =  h  -  X 

Thus, 

du  =  -  dx 

and 

X  =  h  -  u  (A-13d) 

For  the  limits,  note  from  Equation  A-13b  that  when 

x  =  0,  u  =  h  (A-13e) 


(A- 13b) 

(A- 13c) 
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and  when 


f  w, 1  . -  -  ■■'■»!■ '  ■  ■»—•.•*«■■- —  ■- 


T  =  h,  u  =  0 

By  using  Equations  A- 13c  through  £  in  A* 13a  one  obtains 

0 

Sqj  =  -  j ' (a  +  b  h  -  b  u)  e  ^UW  sin  Wj  u  du 

h 

which  can  be  written  as 


(A-13f) 


(A-I4a) 


h 

=  (a  +  b  h)  J' e  ^u)u  sin  uj^  u  du 
0 


sin  ui^  u  du 


(A- 14b) 


These  are  tabulated  integrals,*  and  their  evaluation  is  straightforward, 
but  lengthy.  After  considerable  simplification,  one  obtains 

S01  =  °20  y0  +  a21  yl  (A'l5a) 

where 

"20  5  ir  k  *  +  nX  +  ITh  C(2  ?  '  in  +  2  C  n(X-  l)]j  (A- 15b) 

“21  5  ~k  jn  +  ^Th  1(2  C2  -  D»l>  ♦  2  C  Q(X  -  D]  J  (A-15c) 

It  should  be  noticed  from  Equations  A-12b  and  A-12c  that  the  same  result 

would  be  found  if  the  integration  were  from  y^  to  y^i  thus,  is  giyen 
by 


*See  Reference  16,  pp  99  and  101. 
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°20  yl  +  ®21  y2 


(A-15d) 


Next,  Equations  A-12  are  used  in  the  definition  of  Cni *  Equation  A-2d, 


which  becomes 


u 

CQ1  =  J (a  ♦  b  t)  e  ^  cos  Wj(h-x)  dx 


(A-16a) 


By  using  the  same  change  of  variable  indicated  in  Equations  A-13,  one 
obtains  similarly 


U 

*  (a  t  b  h)  J' e  ^UJU  cos  to^  u 


(A-l6b) 


b  /u  e  ^IJU  cos  u)^  u  du 


which  as  before  contains  integrals  that  are  tabulated  in  Burrington 
(Ref  16).  Again,  after  considerable  simplification,  the  result  can  be 
written  as 


a22  y0  +  °23  yl 


(A- 17a) 


where 


i  N  *  '  C*  +  Fh  12  *  n  ♦  ‘  (2  ?  •  l><*“  l>l}  (A-17b) 
£  |t  -  ITh  12  t  n  *  -  (2  t2  -  1)(X-  l)]j  (A-17c) 
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TWO-EQUATION  COMPUTING  PAIR  FOR  STRAIGHT  LINE  ACCELERATION 


To  permit  verification  of  the  results,  substitute  Equations  A-15a 
and  A-17a  into  Equations  A-4a  and  A-4b  to  9how  that  correct  values  of 
the  constants  were  derived.  This  yields 


!1  *  V.’  °2  *0  *  V°20  *0  *  “21  V 

and 

*1  =  “5  *0  *  Va20  *0  *  °21  *1> 

t  -7(«22  y0  *  «23  y,) 

Defining  new  constants,  these  can  be  written  as 


(A-18a) 


(A-18b) 


Z1  = 

°1  Zo  +  “2 

zo  + 

a29  *0  +  “30  *1 

(A- 18c) 

zl  = 

"4  Zo  +  a5 

zo  + 

a31  *0  +  a32  *1 

(A-18d) 

where 

=  a3  a20 

(A-18e) 

“30 

£  a3a21 

(A-18f) 

“31 

=  "6  °20 

+  a7 

a22 

(A- 18g) 

a32 

5  a6a21 

+  a7 

a23 

(A-18h) 

Using  the  previously  obtained  values  of  the  indicated  a's,  these  are 
found  to  be 


a0Q  ~  -T~  R2  C2  ”  1  +  C  w  h)  ~  +  (2  C  +  w  h)  X  -  2  £ 1  (A-18i) 

*  <u  h  L  n  J 

®30  =  ~T^  [(1  -  2  t2)  A  -  2  t(X-  1)  -  u.  b]  (A-18J) 
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(A-18k) 


“3.  -«♦•«*-*] 

•»  <A-8l) 

These  are  the  equations  used  in  Reference  3,  and  the  constants 
agree,  except  for  two  typographical  Disprints,  which  are  corrected  in 
the  Fortran  program  listing.  One  might  at  first  think  these  expressions 
adequate  if  both  the  relative  velocity  and  displacement  are  desired. 
However,  such  is  not  the  case  because  the  recursive  process  would  then 
involve  more  calculations  and,  thuB,  a  greater  accumulation  of  error. 

The  equations  are  only  given  here  for  completeness. 


NEW  SINGLE  RECURSIVE  EQUATION  FOR  A  STRAIGHT  LINE  ACCELERATION 
APPROXIMATION 

Now  the  results  of  evaluating  the  integrals  SQ1,  S^2,  and  CQ1 ,  with 
the  acceleration  taken  as  a  straight  line  between  the  digitized  values, 
are  obtained  by  merely  substituting  Equations  A-15  and  A-17  into 
Equation  A-9c.  The  result  can  be  arranged  as  follows 

z2  =  a!5  Zo  +  a!6  Z1  +  a24  *0  +  °25  *1  +  a26  *2  (A'19a) 

where  the  new  constants  are  found  to  be 

a24  2  [2  t  X  -  e‘2tu)h(2  C  +  w  h)  +  (1  -  2  £2)  Aj  (A-19b) 

a25  '  ~T^  [“  h  X  ‘  (1  '  1  t2)  I  “  CU  *  e‘2tU,h)]  (A-I9c) 

cr26  s  “T^[(1  -  2  C2)  %  +  2  CU  -X)  -  u>  hj  (A-19d) 
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and  the  first  two  constants  are  repeated  here  for  convenience 


f 


2  I 

%  f 


"l4  ■  2X 


(A-lOh) 

(A-lOi) 


It  is  hoped  that  the  reader  has  observed  that  the  indices  0,  1,  and 
2  on  the  z  and  y  values  do  not  restrict  this  result  to  the  first  three 
points  of  the  data  list,  i.e.,  they  are  perfectly  general;  in  fact,  to 
better  conform  to  other  writings  on  recursive  calculating  equations, 
Equation  A* 19a  can  be  written  for  calculating  the  ith  value  as 

*i  *  °15  *1-2  *  “l4  2i-l  *  °24  >1-2  *  “25  >1-1  +  “24  >1 


Though  it  is  a  trivial  obaervatien,  especially  after  going  through  the 
algebra  three  tines,  the  reader  will  note  on  his  first  time  through  that 


r  * 


2  .  -2twh 


(A-19f) 


E 


is  a  group  that  recurs. 

The  velocity  equation  for  calculating  velocities  from  the  computed 
displacements  is  similarly  obtained  by  substituting  Equations  A-15  and 
A-17  into  A-8b  which  again,  after  simplification,  can  be  wrixten  as 


2i  = 


°11Z  i-l+  °12  zi  +  "27  *i-l  +  °28  h 


(A-20a) 


where 


*27 


28 


2 

ti»  h 


2 

uT  h 


.  n  (2  C  ♦  u>  h)  -  2  t2  ♦  lj 


J  2  t  0  «"2Cwh  ♦  (w  h  -  2  t)  -  C  'I 


(A-20b) 

(A-20c) 


and  from  the  previous  results 


or 


11 


4-  n  w 

<1» 


e-2Cu.h 


(A- 10a  ) 


(A- 10c) 


Equations  A-19e  and  A-20a  are  the  new  result.  Equation  A-I9e  is 
the  equation  that  marches  through  the  data  yielding  relative  displace¬ 
ments,  z^'s,  from  only  previously  computed  z^'s  and  excitation  values. 

Note  also  that  each  step  requires  only  five  multiplications  and  four 
additions.  It  is  theoretically  identical  to  using  both  Equations  A-I8c 
and  A-I6d,  but  avoids  relative  velocity  as  an  intermediate  result. 
Equation  A-20a  is  the  velocity  calculation;  it  exactly  calculates  the 
velocity  at  any  desired  point  from  computed  displacements  and  input 
values . 

Incidentally,  it  is  more  straightforward  to  deduce  Equations  A-l8c 
and  A-18d  directly  from  the  original  differential  Equation  5  by  starting 
with  the  straight  line  acceleration  approximation  in  the  right-hand  side 
and  obtaining  the  needed  particular  solution  by  undetermined  coefficients. 
This  is  the  technique  of  Reference  3.  Then  Equations  A-18c  and  A-18d 
are  transformed  to  Equation  A-19e  by  going  through  the  steps  of 
Equations  A-7  to  A-9  for  this  specific  case  rather  than  for  the  general 
case. 


FURTHER  COMMENTS  ON  ACCURACY  AND  REMAINING  PROBLEMS 

Some  comments  were  made  in  the  report  body  on  accuracy,  and  they 
are  applicable  here.  Additionally,  one  might  say  that  if  both  the 
relative  displacement  and  the  relative  velocity  are  desired,  one  might 
just  as  well  use  the  two-equation  approach  of  References  3  and  10.  That 
would  be  less  accurate  because  the  propagation  of  the  solution  through 
the  data  would  require  eight  as  opposed  to  five  multiplications  per  step 
and  six  as  opposed  to  four  additions  per  step.  Accuracy  would  suffer. 
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One  can  certainly  envision  parabolic  or  cubic  approximation  of  the 
acceleration  input  during  the  integration,  and  cne  would  certainly 
expect  this  to  be  an  improved  approximation.  It  would  permit  one  to 
sample  more  coarsely,  perhaps  at  only  five  points  per  cycle  of  highest 
frequency  present.  But  in  a  sense  this  becomes  defeating  and  introduces 
complications.  In  the  method  presented  here,  the  results  are  only 
evaluated  at  each  point.  As  long  as  the  points  occur  at  least  10  points 
per  cycle,  the  maximum  error  can  only  be  5%.  Making  the  sampling  coarser 
would  then  require  an  additional  procedure  of  evaluating  the  result  in 
between  points.  This  is  suggested  in  Reference  10  and  implemented  in 
Reference  5.  At  this  point  it  appears  to  be  a  substantial  effort  to 
implement  and  test;  the  value  of  the  outcome  is  at  best  uncertain. 

Finally,  what  is  most  disappointing  and  complicating  is  the  fact 
that  all  of  the  constants  in  all  of  the  methods  have  insidious  low 
frequency  problems.  An  examination  of  all  the  a's  reveals  that  crucial 
ones  in  all  of  the  methods  approach  indeterminate  forms  for  very  small 
values  of  tub.  Thus,  one  cannot  without  extreme  caution  evaluate  for 
frequencies  outside  the  range 

10  £  fg/f  £  8,000  (A-21) 

This  range  was  determined  empirically  for  a  computer  that  carries  about 
14  significant  figures.  Double  precision  improves  the  situation,  but 
not  considerably.  Undoubtedly  this  is  being  resolved  by  those  producing 
specialized  shock  spectrum  computing  machinery;  however,  their  algorithms 
are  invariably  proprietary  and,  thus,  not  subject  to  technical  review. 

What  will  have  to  be  done  is  a  careful  series  expansion  of  the  a's,  and 
then  the  programs  will  have  to  shift  to  the  expanded  a's  when  frequencies 
are  low  enough  to  permit  their  use.  Until  such  tiise  as  this  can  be 
done,  the  above  limits  must  be  imposed.  It  is  also  presuaied  that  since 
none  of  the  authors  point  thiB  out,  they  are  unaware  of  the  problem.  It 
is  insidious  because  there  is  no  way  you  can  realize  it  is  happening. 

All  of  a  sudden  spurious  peaks  and  valleys  occur  in  the  spectrum  with  no 
continuity  from  one  frequency  to  the  next,  and  one  finds  that  the  computer 
is  generating  incorrect  values  for  the  a's. 
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LIST 


80/09/23.  19.03.28. 

PROGRAM  SPCTRM 

00100  PROGRAM  SPCTRMI INPUT, OUTPUT. TAPE! .TAPE10) 

00110  DIMENSION  YI2050 > , Z ( 2050) , 0UTI5. 500) 

00120  PRINT, *INPUT  ACCELS  IN  G'S* 

00130  PRINT, *NO. OF  SAMPLES  OF  INPUT  DATA  (I)* 

00140  READ,N1 

00150  PRINT, ★APPROX  LON  FREQ,  (HZ),  (F)* 

00160  READ, FLOW 

00170  PRINT, ★APPROX  HIGH  FREQ,  <F>* 

00180  READ,FHIGH 

00190  PRINT, ★SAMPLING  RATE  (SAMPLES/SEC.)  (F)* 

00200  READ,S 

00210  PRINT, ★FREQS  PER  DECADE,  (AST  25),  (F)* 

00220  READ.FPD 

00230  PRINT, ★SCALING  FACTOR  (G,S/DATA  UNIT)  (F)* 

00240  READ,SF 

00250  PR  I NT, ★ZERO  OFFSET  (DATA  UNITS)  (F)* 

00260  READ,ZO 

00270  FLOWLOG-ALOGIO(FLOW) 

00280  C2-FL()AT(  INT(FL0WL0G)-1  ) 

00290  AJL()W«FPD*(FL0WL0G-C2) 

00300  JSTART® I NT ( A  JLOW ) 

00310  I F ( A JLOW . NE . FLOAT ( JSTART ) ) JSTART- JSTART* 1 
00320  FSTART- 1 0. ★*( FLOAT < JSTART ) /FPD+C2 ) 

00330  JSTOP- I NT ( FPD* ( ALOG 1 0 ( FH IGH) -C2 ) ) 

00340  FSTOP- 1 0 . ★★ ( FLOAT <  JSTOP ) /FPD+C2 ) 

00350  NFREQS-JST0P-JSTART+1 

00360  PRINT, ★LOWEST  AND  HIGHEST  FREQUENCIES  ARE*, FSTART, FSTOP 

00370  PRINT, *N().  OF  FREQS  IS*,NFREQS 

00380  TP  I -6. 283 1 853 1 

00390  READ (1,69)(Y(I),I=I,NI ) 

00400  69  F0RMAT(5X,4( IX, 15)) 

00410  GSF«SF*386.088 

00420  DO  33  1*1 ,N I 

00430  Y ( I ) -GSF* ( (FLOAT ( Y ( I  ) ) ) -Zo ) 

00440  33  CONTINUE 

00450  PRINT, *DAMPING  RATIO?  (F)* 

00460  READ,ZETA 
00470  H*1 ,/S 

00480  ETA*SQRT< I . -ZE TA*ZETA ) 

00490  GI*2.*ZETA 
00500  G2»l .-GI*ZETA 
00510  DO  111  J-JSTART, JSTOP 
00520  F»IO. ★★(FLOAT! J)/FPD*C2) 

00530  WOM-TPI*F 
00540  G3*W0M*H 
00550  G4-EXP(-ZETA*G3) 

00560  A I 5*-G4*G4 
00570  G5-ETA*G3 
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00580  G6-G4/ETA*SIN(05) 

00590  CHI«G4*CoS(G5> 

00600  G7-CHI/G6 
00610  G8—AI5/G6 
00620  A 1 6-2. *CH I 
00630  A I  l*-WOM*G8 
00640  AI2»W()M*(G7-ZETA) 

00650  G9-G3*W<)M 
00660  010*01 *G3 
00670  0 1 t "02*06 
00680  0I2*G9*N()M 

00690  A24*(GI*CHI4A15*GI0+G1 I )/GI2 

00700  A25-2./GI2*(G3*CHI-GII-ZETA*(1.+A!5)) 

00710  A26- (G 1 1  +G I *(  l .-CHI )-G3)/G!2 
00720  A27-(G1*G7-G8*010*G2)/G9 
00730  A28-(01*GB-MG3-G1  )*(G7-ZETA)-1 .  )/G9 
00740  ZIMIN-0. 

00750  ZI MAX-0. 

00760  Z( I )-A26*Y ( I ) 

00770  Z(2)-AI6*Z( I )*A25*Y< ! )*A26*Y(2) 

00780  DO  1 3  1-3, N I 

00790  Z ( I ) *AI 5*Z ( I -2 ) +A 1 6*Z ( I- 1 ) +A24*Y ( 1-2 )♦ A25*Y ( I- 1 ) *A26*Y ( I ) 
00800  IF(Z( I ) .GT.Z IMAX )Z I  MAX-Z ( I ) 

00810  IF(Z(I).LT.ZIMIN)ZIMIN-ZCI) 

00820  13  CONTINUE 
00830  ZO-Z(NI). 

00840  ZDO-AI I*Z(N1 -I )+Al 2*Z(NH *A27*Y CN I - 1 ) +A28*Y(N1 ) 

00850  CALL  RESID(WOM, ZETA.ETA.ZO, ZDO.ZRMAX, ZRMIN > 

00860  ZMIN-AMINI (ZIMIN, ZRMIN) 

00870  ZMAX-AMAX I (Z I  MAX , ZRMAX ) 

00880  SV MIN-W()M*ZM  IN 
00890  SVMAX«WOM*ZMAX 
00900  OUT(l ,J)-Z MIN 
00910  ()UT(2,J)-ZMAX 
00920  ()UT<3, J)-SVMIN 
00930  ()UT(4, J)-SVMAX 
00940  ()UT(5, J)*F 
00950  IJ1  CONTINUE 
00960  DO  410  I-JSTART, JSTOP 

00970  410  WRITE! 10,425  >0UT<5, 1 ) , (OUT ( KK, I ) , KK- I ,4 ) 

00980  PRINT,/*  PLEASE  NOTE*/ 

00990  PRINT, *0N  THE  OUTPUT  FILE  THE  DATA  COLUWS  READ  AS* 

01000  PRINT. *FREOUENCY( HZ),  ZMIN(IN),  ZMAX,  SVMIN(IPS),  SVMAX* 

01010  PRINT, *TYPE  "SAVE, TAPE  1 0- _ "* 

01020  PRINT, *TYPE" RETURN, TAPE  10"* 

01030  425  F()RMAT(F  12.5, 4E 15. 6) 

01040  END  FILE  10 
01050  STOP 
01060  END 

01070  SUBROUTINE  RES ID(WOM, ZETA, ETA, Zo, ZDO.ZRMAX, ZRMIN ) 

01080C  THIS  SUBROUTINE  COMPUTES  RESIDUAL  RESPONSE 
01090  B-ZO 

01100  DELT-ASIN(ZETA) 

OHIO  A-(ZO*ZETA*ZDO/WOM)/ETA 
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Oi 120  PI-3.141592654 
01130  IF(A.EQ.O. .AND.B.EQ.O. ) Go  TO  24 
01140  IF(A.E0.0.)  00  TO  26 
01150  IFIB.EQ.O. )  00  TO  28 

01160  IF (A. AND.B.OT.O.  .OR . A. AND. B . LT.O. )  GOTO  30 
01170  BETA 1 »P I -ATAN < -A/B ) 

01180  00  TO  80 

01190  30  BETA I -ATAN (A/B) 

01200  GO  TO  80 
01210  28  BETA! -PI/2. 

01220  GO  TO  80 
01230  24  ZRMIN-O. 

01240  ZRMAX-O. 

01250  00  TO  100 
01260  26  BETA 1-0. 

01270  80  IF(BETA1.LT.DELT)BETA1-BETAI+PI 
01280  NOT I -BETA 1 -DELT 
01290  HDT2-MDTH-PI 

01300  Z1-EXP(-ZETA*HDT1/CTA)*(A*SIN<WDT1 )*B*C()S(WDTl >) 
01310  Z2-EXP (-ZETA*HDT2/ETA )*<A*$IN< NDT2 ) +B*COS ( HDT2 ) ) 
01320  ZRMAX-AMAXI (Z1.Z2) 

01330  ZRMIN-AMINI ( Z I tZ2) 

01340  100  RETURN 
QI350  END 


80/09/23.  19.06.03. 

PROGRAM  PLTVELF 

00100  PROGRAM  PLTVELF( INPUT. OUTPUT, TAPEI ,TAPE1 0) 

00110  DIMENSION  FREQ<250) , VEL(250) , IFRE0(250) , IVELC250) , VELN<250) 

00120  PRINT,*TYPE#FTV#FOR  A  FT  PLOT  OR  TYPE  #SS#  FOR  A  SS  PLOT* 

00130  READ73.FT 
00140  73  FORMAT (A3) 

00150  IF(FT.E0.3HFTV)G0TO74 

00160  PRINT,*TYPE#*l#FOR  P()S.SPECT.OR#-l#FOR  NEG.SPECT.OR#*0#FOR  O.ASS* 
00170  READ  99, IPLOT 
00180  99  FOR MAT ( 12 ) 

00190  74  PRINT, *N0. OF  AMPLITUDES  TO  BE  PLOTTED* 

00200  READ.N 

00210  PRINT, *TYPE#PLTL#F()R  A  LINE  PLOT  ORfPLTP#  FOR  A  PT.PLOT* 

00220  READ  7 I, PLOT 
00230  71  F()RMAT(A5) 

00240  IF (FT.E0.3HFTV)GoT()75 
00250  IF ( IPL0T>4, 5, 6 

00260  75  READ<»,20HFREQ<J),VEL<J),J«I  ,N) 

00270  20  FORMAT(2E15.6) 

00280  GO  TO  7 

00290  4  READ! I , 10) (FREQ< J) ,VEL(J), J«l ,N> 

00300  18  FORMAT(F12.5,30X,EI5.6) 

00310  GO  TO  7 

00320  17  F0RMAT(FI2.5,45X.EI5.6) 

00330  6  READ(l,l7)(FRE0(J),VELiJ),J-l,N) 

00340  7  K-0 

00350  IF(FREO(l).LT.I.)  K-2500 

00360  I F ( K . EQ . 0 ) PR  I NT, *THE  ORIGIN  IS  AT  A  FREQUENCY  OF  KHZ)* 

00370  IF(K .NE.OPRINT, *THE  ORIGIN  IS  AT  A  FREOUENCY  OF  0.1  HZ.* 

00380  DO  88  I-I.N 

00390  IFREQ(  I  MNT<ALOGI0<FREQ<  I  ))*9999./4.) 

00400  I FREQ ( I )■ IFREQ( I )*K 
00410  VEL( I )-ABS(VEL( I ) ) 

00420  88  I VEL( I )-INT( (  ALOGI 0(VEL< I >  >  +  1 . )*9999./5. ) 

00430  PRINT, *ALL  THE  VALUES  ARE  IN* 

00440  PRINT, *SET  UP  PLOTTER  TO  PLoTl  TURN  ON  PUNCHl  HIT  CR* 

00450  PAUSE 

00460  PRINT, *HERE  IS  SOME  LEADER* 

00470  PRINT7I  ,PL()T 

00480  DO  30  1-1, N 

00490  30  PRINT  40, I FREQ < I > , IVEL( I ) 

00500  40  FORMATI2I5) 

00510  PRINT,*PLTT* 

00520  PRINT. *SHUT  OFF  PLOTTER!  TURN  OFF  PUNCH* 

00530  GO  TO  10 

00540  19  FORMATiFI 2.5, 30X.2EI5.6) 

00550  5  READ! 1,19) (FREQ( J) , VELNI J) , VEL< J) , J-1 ,N) 

00560  K-0 

00570  IFIFREQd  ).LT.  I.)  K-2500 

00580  IF(K.EQ.O)PRINT,*THE  ORIGIN  IS  AT  A  FREQ. OF  1.0  <HZ)* 

00590  IF (K.EQ.2500)PRINT,*THE  ORIGIN  IS  AT  A  FREQ. OF  0.1  (HZ)* 

00600  DO  21  I-l ,N 

00610  IFRE0< I )-INT(ALOGI0<FREQ< I ) )*9999./4. ) 
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00620  IFREQ( I )» I FREQ ( I )+K 
00630  VELN<I)»ABS(VELNU)) 

00640  VEL( I )»ABS (VEL( I ) ) 

00650  VEUD-AMAXI  <VELN( I > , VEL( I >) 

00660  21  IVEL<I>*INT<<ALOGIO(VEL(I)>+l.>*9999./5.) 

00670  PRINT, *ALL  THE  VALUES  ARE  IN* 

00680  PRINT, *SET  UP  PLOTTER  TO  PLOT*  TURN  ON  PUNCH*  HIT  CR* 
00690  PAUSE 

00700  PRINT, *HERB  IS  SOME  LEADER* 

00710  PRINT7I  ,  PL()T 
00720  DO  31  I«l ,N 

00730  31  PRINT  40, IFREQI I ) , I VEL< I ) 

00740  PRINT, *PLTT* 

00750  PRINT, *SHUT  OFF  PLOTTER*  TURN  OFF  PUNCH* 

00760  10  STOP 
00770  END 


Appendix  C 


DETAILS  OF  THE  RESIDUAL  SPECTRUM  CALCULATION 


The  table  given  in  the  nain  text  ia  repeated  below  to  organize  the 
exaaination  of  the  nine  cases 


A  >  0 

A  =  0 

A  <  0 

B  >  0 

I 

II 

III 

B  =  0 

IV 

V 

VI 

B  <  0 

VII 

VIII 

IX 

Because  many  values  of  0  will  satisfy  Equation  15a,  one  iaagines  the 
signs  of  A  and  B  in  Equation  15  to  deduce  the  proper  value  of  S  for  0j. 


Case  I,  A  >  0.  B  >  0 

Consider  a  sketch  of  A  cos  0  and  B  sin  0  as  shown  in  Figure  C-l. 


As  can  be  seen,  is  an  acute  angle,  with  tangent  defined  by  Equation  15a, 
and  using  principal  angle  notation 
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(C-l) 


Pj  =  Tan”1  (A/B) 

Recall  that  (w^t)  must  be  greater  than  zero;  thus,  from  Equation  14,  the 
following  check  is  made 

Pj  -  6  £  0  (02) 

If  Condition  02  is  not  fulfilled,  the  value  of  p^  must  be  increased  by 
n.  The  value  of  p^  is  given  by  Equation  16  in  all  cases. 

Case  II.  A  •  0,  B  >  0 

Consider  Figure  02. 


Figure  02. 

In  every  case 

Px  *  0  (03) 

However,  in  every  case  with  damping,  this  value  of  Pj  will  not  satisfy 
Condition  02;  thus, 

Pj  =  n,  if  C  =  0  (04) 

Again  Pj  will  be  given  by  Equation  16. 
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Consider  the  sketch  of  Figure  C-3. 


Clearly  0j  is  betweea  n/2  aud  n,  but  note  that  if  the  negative  of  A  cos 
0  is  drawn  dotted,  an  angle  0^  can  be  defined  as 

=  Tan'1  (-A/B)  (C-5) 

The  angle  0j  it  then  given  by  (n  -  0p  or 

0j  =  n  -  Tan'1  (-A/B)  (C-6) 

Condition  C-2  would  again  be  used  to  see  if  0j  should  be  increased  by  n 
and  Equation  16  would  give 

Case  IV.  A  >  0.  B  =  0 

Consider  the  sketch  of  Figure  C-4. 


Figure  C-4. 


Clearly  pj  =  n/2.  Condition  C-2  would  be  checked,  and  Equation  16  used 
for  p2. 

Case  V.  A  =  0,  B  =  0 

Going  back  to  the  definition  of  A  and  B  in  Equations  11a  and  lib, 

one  sees  that  both  tQ  and  zQ  are  zero.  Hence,  no  residual  response 

results,  and  both  z  .  and  z  are  zero. 

min  max 

Case  VI,  A  <  0,  B  =  0 

Consider  the  sketch  of  Figure  C-5. 


Figure  C-5. 

Clearly  =  n/2.  Condition  C-2  would  be  checked,  and  Equation  16  used 
for  P2. 

Case  VII,  B  <  0,  A  >  0 

Consider  the  sketch  of  Figure  C-6. 
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In  this  case  pj  lies  between  n/2  and  n.  If  the  negative  of  B  is  used, 
the  acute  angle  p'  can  be  obtained  and  subtracted  from  n  as  indicated 
below. 

P'  =  Tau*1  (A/-B) 
p,  =  n  -  p’ 

Case  VIII.  A  =  0.  B  <  0 

Consider  this  situation  drawn  in  Figure  C-7. 


As  in  Case  II, 

Pj  =  0  (C-8a) 

but  again  in  every  case  with  damping,  this  value  of  pj  will  not  satisfy 
Condition  C-2;  thus, 

Px  =  n,  if  C  =  0  (C-8b) 

As  in  all,  cases  P2  is  given  by  Equation  16. 


(C-7a) 

(C-7b) 
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Figure  C-8. 

This  use  is  the  inverse  of  Case  I,  and  pj  will  always  be  acute  and  given 
simply  by 

pj  =  Tan-1  (A/B)  (C-9) 

If  Condition  C-2  is  not  fulfilled,  the  value  of  Pj  from  Equation  C-9 
must  be  increased  by  n;  the  value  of  p^  is  obtained  from  Equation  16. 

That  completes  the  examination  of  the  nine  possible  cases  identified 
in  Table  1.  By  going  through  them,  one  can  organize  the  results  as 
follows . 

Case  V:  A  =  B  =  0, 

no  residual  response  occurs 
Cases  I,  IX;  A>0,  B  >  0;  or  A  <  0,  B  <  0, 

Pj  =  Tan'1  (A/B) 

Cases  III,  VII:  A  <  0,  B  >  0;  or  B  <  0,  A  >  0, 

Pj  =  n  -  Tan'1  (-A/B) 

Cases  II,  VIII:  A  =  0,  B  >  0;  or  A  =  0,  B  <  0, 

Pi  =  0 
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Case*  IV,  VI:  A  >  0,  B  =  0;  or  A  <  0,  B  =  0, 


Mow  this  summary  can  be  organized  into  a  more  concise  set  of  rules  for 
Pj  determination  as  follows: 

1.  If  both  A  and  B  are  zero,  no  residual  response  results. 

2.  If  not,  and  A  =  0,  pj  =  0. 

3.  If  not,  and  B  =  0,  Pj  =  n/2. 

4.  If  not,  and  A  and  B  have  the  same  sign, 

pj  =  Tan"1  (A/B) 

5 .  If  not , 

pj  s  n  -  Tan'1  (-A/B) 

Now  given  a  value  of  p^,  one  applies  Condition  C-2,  or  if 
Px  -  6  <  0 

Pi  =  “Void*" 

which  yields  acceptable  values  of  p^  from  which  (u^t)^  can  be  obtained 
with  Equation  14  as  follows 

(u»dt)i  =  Pj  -  6  (C-10) 

The  second  value  of  (lOjt)  is  obtained  from  Equation  16  and  is 
(u*dt)2  =  ^d0!  +  71 
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Appendix  D 


FREQUENCY  GENERATION  FOR  EQUALLY  SPACED  VALUES 
WHEN  PLOTTED  ON  LOGARITHMIC  PAPER 


For  a  nuaber  of  reasons  it  is  rdvantageous  to  plot  shock  spectra  on 
four-coordinate  paper,  which  is  a  logarithmic  paper,  as  discussed  in 
Reference  2.  To  do  this  requires  the  frequencies  to  be  selected  in  the 
following  Banner.  For  equally  spaced  values  of  frequency  on  logarithmic 
paper  one  wants  the  logarithai  of  the  frequency  to  have  equally  spaced 
values.  One  can  obcerve  that,  if  J  is  the  computer  "do-loop"  index, 
that  increaents  by  unity  for  each  step,  the  logs  will  be  equally  spaced 
if 


l°g10F  =  V  +  C2  (D-l) 

Corresponding  to  this,  the  frequency  will  be  given  by 
(C.J  ♦  C.) 

F  =  10.  4  (D-2) 

The  user  will  select  a  lowest  frequency  of  interest,  F^ow,  a  highest 

frequency  of  interest,  F^^,  and  the  nuaber  of  frequencies  per  decade 

at  which  he  wants  shock  spectrum  values  calculated,  F  ..  The  remainder 

pd 

of  this  argument  is  more  simply  written  in  FORTRAN  because  integer  logic 
is  necessary.  Thus,  Equations  D-l  and  D-2  become 

ALOGIO(F)  =  C1*J  +  C2  (D-3) 

F  =  10.**(C1*J  +  C2)  (D-4) 
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Each  time  the  log  of  the  frequency  increases  by  unity,  the  frequency 

will  have  increased  by  a  factor  of  10,  or  gone  through  a  decade;  thus, 

Cl  is  the  reciprocal  of  F  j,  and  Equations  D-3  and  D-4  may  be  written 

A10G  10(F)  =  J/FPi)  +  C2  (D-5) 

F  =  10.**(J/FPD  ♦  C2)  (D-6) 

Nov,  rather  than  start  and  stop  at  the  precise  low  and  high  frequencies 
selected  by  the  user,  one  will  start  at  FSTARI  and  stop  at  FSTOP,  defined 
~to  be  the  frequencies  computed  by  Equation  D-6  from  the  values  of  J 
given  as  JSTART  and  JSTOP.  Using  functional  notation  this  can  be  indi¬ 
cated  by  saying 

FSTOP  =  F (JSTOP) 

FSTART  =  F (JSTART) 

The  variables  JSTOP  and  JSTART  are  defined  such  that 

F (JSTART  -  1)  .LT.  FLOW  .LE.  F( JSTART)  (D-9) 

and 

F( JSTOP)  .LE.  FHIGH  .LT.  F( JSTOP  ♦  1)  (D-10) 

This  is  convenient  because  the  user  will  in  general  select  an  FPD,  an 
FLOW,  and  an  FHIGH  that  are  not  mutually  consistent.  In  this  way,  if 
the  user  selects  integer  values  of  FPD,  it  is  convenient  to  arrange  C2 
so  that  frequency  values  divisible  by  10  are  always  included  in  the  list 
selected. 

In  Equation  D-S,  C2  is  an  integer  used  so  thct  J  does  not  have  to 
assume  zero  or  negative  values.  Whenever  the  starting  frequency  is  less 
than  unity,  the  log  of  F  will  be  negative.  If  C2  is  one  less  than  the 
characteristic  of  the  log  of  FSTART,  J  will  always  have  to  be  greater 
than  zero.  Thus,  define: 


(D-7) 

(D-8) 
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C2 


INT(ALOG10(FLOW))  -  1 


(D-ll) 


To  formulate  a  procedure  to  fiad  JSTART,  proceed  as  follows. 
Substitute  Equation  D-6  into  Condition  D-9;  take  the  log  of  all  three 
terms,  which  yields 

(JSTART  -  1)/FPD  ♦  C2  .LT.  JLOW/FPD  +  C2  .LE. 

(JSTART/FPD  +  C2)  (D-12a) 

Subtract  C2  from  each  group,  and  then  multiply  through  by  the  positive 
number  FPD  which  will  yield 

JSTART  -  1  .LT.  JLOW  .LE.  JSTART  (D-12b) 

where  JSTART  is  an  integer  and  JLOW  is  floating  point  and  formed  from 
Equation  D-5  with  F  equal  to  FLOW  or 

JLOW  =  FPD*(ALOG10(FLOW)  -  C2)  (D-13) 

In  general,  JLOW  will  not  have  integral  value;  for  the  caae  where  JLOW 
is  not  integral 

JSTART  =  INT(JLOW)  +  1  (D-14) 

but  when  JLOW  does  have  integral  value 

JSTART  =  I NT (JLOW)  (D-15) 

This  can  be  programed  as  follows 

FLOWLOG  «  ALOGIO(FLOW) 

C2  =  FLOAT (I NT (FLOWLOG)  -  1) 

AJLOW  =  FPD* (FLOWLOG  -  C2) 
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JSTART  =  IHT(AJIOW) 

IFCAJLOW  .ME.  FLOAT (J8TART))  JSTART  =  JSTART  +  1 
FSTART  =  10.**(FLOAT(JSTART)/FPD  ♦  C2) 

A  few  values  are  hand-coaputed  to  docuaent  that  this  does  deliver 
Condition  D-9.  The  results  are  shown  in  Table  D-l. 

To  derive  formulas  for  JSTOP  and  FSTOP ,  substitute  Equation  D-6  in 
Condition  D-10  and  take  the  log  of  all  three  tents,  which  yields 

JSTOP/FPD  ♦  C2  .LE.  JHIGK/FPD  ♦  C2  .LT.  (JSTOP  +  1)/FPD  +  C2 

Subtracting  C2  from  each  tera,  and  then  aultip lying  each  by  the  positive 
number  FPD,  yields 

JSTOP  .LE.  JHIGH  .LT.  JSTOP  +  1  (D-16) 

where  JSTOP  is  an  integer,  and  JHIGH  is  found  froa  Equation  D-5  with  F 
equal  to  FH1GH,  or 

JHIGH  =  FPDKALOGIO(FHIGH)  -  C2)  (D-17) 

In  general,  JHIGH  will  not  be  integral,  but  even  when  it  is,  Condition 
D-16  will  be  satisfied  if  we  tske 

JSTOP  =  INT (JHIGH)  (D-18) 

Therefore,  to  get  JSTOP  and  FSTOP,  one  adds  the  following  steps  to  the 
prograa  given  after  Equation  D-1S 

JSTOP  *  IHT(FPD*(ALOG10(FHIGH)  -  C2)) 

FSTOP  *  10.**(FLOAT(JSTOP)/FPD  ♦  C2)  (D-19) 
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CA;  PWO  Newport  RI,  PWO  Portsmouth.  VA;  SCE  (D.  Kaye);  SCE  San  Diego,  CA;  SCE,  Camp 
Pendleton  CA;  SCE,  Guam;  SCE,  Oakland  CA 
NAVSCOLCECOFF  C35  Port  Hueneme,  CA;  CO,  Code  C44A  Port  Hueneme.  CA 
NAVSEASYSCOM  Code  05M13  (Newhouse)  Wash  DC;  Code  62T2.  Wash  DC;  SEA-04H,  Wash  DC 
NAVSEC  Code  6034  (Library),  Washington  DC 

NAVSECGRUACT  PWO,  Adak  AK;  PWO,  Edzell  Scotland;  PWO,  Puerto  Rico;  PWO,  Torn  Sta,  Okinawa 
NAVSHIPREPFAC  SCE  Subic  Bay 

NAVSHIPYD;  Code  202.4,  Long  Beach  CA;  Code  202.5  (Library)  Puget  Sound,  Bremerton  WA;  Code  380, 
(Woodroff)  Norfolk.  Portsmouth.  VA;  Code  400,  Puget  Sound;  Code  400.03  Long  Beach,  CA;  Code  404 
(LT  J.  Riccio),  Norfolk.  Portsmouth  VA;  Code  410.  Mare  Is.,  Vailcjo  CA;  Code  440  Portsmouth  NH;  Code 
440.  Norfolk;  Code  440.  Puget  Sound,  Bremerton  WA;  Code  450,  Charleston  SC:  Code  453  (Util  Supr), 
Vallejo  CA;  L.D.  Vivian;  Library.  Portsmouth  NH;  PWD  (Code  400),  Philadelphia  PA;  PWO.  Mare  Is.; 
PWO,  Puget  Sound:  SCE.  Pearl  Harbor  HI:  Tech  Library.  Vallejo,  CA 
NAVSTA  CO  Naval  Station.  Mayport  FL;  CO  Roosevelt  Roads  P  R.  Puerto  Rico;  Dir  Mech  Engr.  Gtmo; 

Engr.  Dir  .  Rota  Spain;  Long  Beach.  CA;  Maim.  Corn.  Div..  Guantanamo  Bay  Cuba,  Maim.  Div.  Dir/Codc 
531,  Rodman  Canal  Zone;  PWD  (LTJG.P.M.  Motolenich),  Puerto  Rico;  PWO  Midway  Island.  PWO, 

Keflavik  Iceland;  PWO.  Mayport  FL;  ROICC.  Rota  Spain;  SCE,  Guam;  SCE.  San  Diego  CA;  SCE,  Subic 
Bay,  R.P.;  Utilities  Engr  Off.  (A.S.  Ritchie),  Rota  Spain 
NAVSUBASE  ENS  S.  Dove.  Groton.  CT;  SCE.  Pearl  Harbor  HI 

NAVSUPPACT  CO.  Seattle  WA;  Code  *. ,  12  Marine  Corps  Dist.  Treasure  Is..  San  Francisco  CA;  Code  413. 

Seattle  WA:  LTJG  McGarrah,  SEC.  Vallejo,  CA;  Plan/Engr  Div.,  Naples  Italy 
NAVSURFWPNCEN  PWO.  White  Oak,  Silver  Spring.  MD 
NAVTECHTRACEN  SCE.  Pensacola  FL 

NAVWPNCEN  Code  2636  (W.  Bonner),  China  Lake  CA;  Code  3276  China  Lake.  CA;  PWO  (Code  26).  China 
Lake  CA.  ROICC  (Code  702),  China  Lake  CA 

NAVWPNSTA  (Clebak)  Colts  Neck,  NJ.  Code  09.  Concord.  CA:  Code  092,  Colls  Neck  NJ;  Code  092A  (C. 
Fredericks)  Seal  Beach  CA;  Maim.  Control  Dir  ,  Yorktown  VA;  Naval  Weapons  Handling  Center.  Earle, 

NJ 

NAVWPNSTA  PW  Office  (Code  09CI)  Yorktown.  VA 
NAVWPNSTA  Safety  Br,  Concord.  CA;  Security  Offr,  Colts  Neck  NJ 
NAVWPNSUPPCEN  Code  09  Crane  IN 
NCBU  405  OIC,  San  Diego,  CA 

NCBC  Code  10  Davisville,  RI;  Code  155,  Port  Hueneme  CA;  Code  156,  Port  Hueneme.  CA:  Code  400. 

Gulfport  MS;  PW  Engrg.  Gulfport  MS,  PWO  (Code  8(1)  Port  Hueneme.  CA;  PWO,  Davisville  RI 
NCBU  411  OIC,  Norfolk  VA 
NCR  20,  Commander 
NCSO  BAHRAIN  Security  Olfr.  Bahrain 

NMCB  5.  Operations  Dept.;  74,  CO;  Forty,  CO;  THREE,  Operations  Off 

NORDA  Code  440  (Ocet  Rsch  Off)  Bay  St  Louis  MS 

NRL  Code  8400  Washington.  DC;  Code  8441  (R  A.  Skop).  Washington  DC 

NSC  Code  54.1  (Wynne).  Norfolk  VA 

NSD  SCE.  Subic  Bay,  R.P. 

NTC  Commander  Orlando.  FL 
NTSB  Chairman,  Wash  DC 

NUSC  Code  131  New  London,  CT;  Code  EAI23  (R  S  Munn),  New  London  Cl 
OCEANSYSLANT  LT  A.R.  Giancola.  Norfolk  VA 

OFFICE  SECRETARY  OF  DEFENSE  OASD  (MRA&I.)  Pentagon  (T  Cjshcrjj).  Washington.  DC 
ONR  Code  700F  Arlington  VA 
PHIBCB  I  P&E.  Coronado.  CA 
PMTC  Pat.  Counsel.  Point  Mugu  CA 

PWC  (l.t  E  S.  Agonoy)  Pensacola,  FL;  ACE  Office  (LTJG  St.  Germain)  Norfolk  VA.  CO  Norfolk,  VA.  CO. 
(Code  10).  Oakland.  CA:  CT).  Great  Lakes  II.;  Code  10.  Great  lakes.  I!.,  Code  I2ti,  Oakland  C'A,  Codr 
1 21X7 .  (Library)  San  Diego,  CA;  Code  128,  Guam,  Code  154,  Great  Lakes.  It..  Code  21*1.  Great  Lakes  IL, 
Code  220  Oakland.  CA;  Code  220.1.  Norfolk  VA;  Code  30C.  San  Diego.  C'A.  Code  44*1,  Great  Lakes.  IL; 
Code  400,  Oakland.  CA;  Code  44*1.  Pearl  Harbor.  HI;  Code  44*),  San  Diego.  CA;  Code  420,  Great  Lakes. 
IL;  Code  420,  Oakland.  CA;  Code  42B  (R  Pascua).  Pearl  Harbor  HI,  Code  54>5A  (II  Wheeler);  Code  64*), 
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Great  Lakes.  IL;  Code  601.  Oakland.  CA;  Code  610,  San  Diego  Ca;  Code  700.  Gteat  Lakes,  1L;  Code  700. 
San  Diego.  CA.  LTJG  J.L.  McClatne.  Yokosuka,  Japan;  Utilities  Officer.  Guam;  XO  (Code  20)  Oakland, 

CA 

SPCC  PWO  (Code  120)  Mechanicsburg  PA 
NAF  PWO  (Code  30)  El  Centro.  CA 

U  S.  MERCHANT  MARINE  ACADEMY  Kings  Point.  NY  (Reprint  Custodian) 

USCG  (Smith).  Washington,  DC;  G-EOE-4/61  (T.  Dowd),  Washington  DC 
USEUCOM  (ECJ4T-LO).  Wright,  Stuttgart.  GE 

USNA  Ch.  Mech  Engr.  Dept  Annapolis  MD;  PWD  Engr.  Div.  (C.  Bradford)  Annapolis  MD 

California  state  university  long  beach,  ca  (CHelapatd 

CORNELL  UNIVERSITY  Ithaca  NY  (Serials  Dept,  Engr  Lib.) 

DAMES  &  MOORE  LIBRARY  LOS  ANGELES.  CA 
ILLINOIS  STATE  GEO  SURVEY  Urbana  1L 

LEHIGH  UNIVERSITY  Bethlehem  PA  (Fritz  Engr.  Lab  No.  13.  Beedle);  Bethlehem  PA  (Linderman  Lib. 

No. 30.  Flecksteiner) 

MICHIGAN  TECHNOLOGICAL  UNIVERSITY  Houghton,  Mi  (Haas) 

MIT  Cambridge  MA;  Cambridge  MA  (Rm  10-500,  Tech.  Reports.  Engr  Lib  );  Cambridge  MA  (Whitman) 

NEW  MEXICO  SOLAR  ENERGY  INST.  Dr  Zwibel  Las  Cruces  NM 
NY  CITY  COMMUNITY  COLLEGE  BROOKLYN.  NY  (LIBRARY) 

PURDUE  UNIVERSITY  Lafayette,  IN  (CE  Engr  Lib) 

CONNECTICUT  Hartford  CT  (Dept  of  Plan.  &  Energy  Policy) 

SEATTLE  U  Prof  Schwaegler  Seattle  WA 

SOUTHWEST  RSCH  INST  R  DeHart,  San  Antonio  TX 

STANFORD  UNIVERSITY  Engr  Lib.  Stanford  CA 

STATE  UNIV  OF  NEW  YORK  Buffalo.  NY 

TEXAS  A&M  UNIVERSITY  W  B  Ledbetter  College  Station.  TX 

UNIVERSITY  OF  CALIFORNIA  BERKELEY.  CA  (CE  DEPT.  GERWICK);  Berkeley  CA  (E  Pearson), 
DAVIS,  CA  (CE  DEPT.  TAYLOR);  LIVERMORE.  CA  (LAWRENCE  LIVERMORE  LAB.  TOKARZ) 
UNIVERSITY  OF  DELAWARE  Newark.  DE  (Dept  of  Civil  Engineering.  Chesson) 

UNIVERSITY  OF  HAWAII  Honolulu  HI  (Dr  Szilard) 

UNIVERSITY  OF  ILLINOIS  Metz  Ref  Rm.  Urbana  IL.  URBANA.  IL  (LIBRARY);  URBANA.  IL 
(NEWMARK);  Urbana  IL  (CE  Dept,  W.  Gamble) 

UNIVERSITY  OF  MASSACHUSETTS  (Heroncmus),  Amherst  MA  CE  Dept 
UNIVERSITY  OF  MICHIGAN  Ann  Arbor  Ml  (Richart) 

UNIVERSITY  OF  NEBRASKA-LINCOLN  Lincoln.  NE  (Ross  Ice  Shelf  Proj  ) 

UNIVERSITY  OF  NOTRE  DAME  Katona.  Notre  Dame.  IN 

UNIVERSITY  OF  TEXAS  Inst.  Marine  Sci  (Library),  Port  Arkansas  TX 

UNIVERSITY  OF  TEXAS  AT  AUSTIN  AUSTIN  TX  (THOMPSON);  Austin.  TX  (Breen) 

UNIVERSITY  OF  WASHINGTON  Dept  of  Civil  Engr  (Dr  Mattock),  Seattle  WA;  SEATTLE.  WA 
(MERCHANT) 

UNIVERSITY  OF  WISCONSIN  Milwaukee  Wl  (Cir  of  Great  Lakes  Studies) 

AMMAN  &  WHITNEY  CONSULT  ENGRS  N  Dobbs  New  York,  NY 
ARV1D  GRANT  OLYMPIA.  WA 
ASSOC  AMER  RR  Bureau  of  Explosive,  Chicago.  IL 
ATLANTIC  RICHFIELD  CO  DALLAS,  TX  (SMITH) 

AUSTRALIA  Dept  of  Const,  Canberra;  Dept  of  Hsng  &  Const,  Sydney 
BECHTEL  CORP  SAN  FRANCISCO.  CA  (PHELPS) 

BROWN  &  CALDWELL  E  M  Saunders  Walnut  Creek,  CA 

CANADA  Mem  Univ  Newfoundland  (Chari),  St  Johns,  Trans-Mnt  Oil  Pipe  Lone  Corp.  Vancouver,  BC  Canada 
COLUMBIA  GULF  TRANSMISSION  CO  HOUSTON.  TX  (ENG  LIB  ) 

HERCULES  INC  D  Richardson,  Magna.  UT 
DURLACH.  O  NEAL.  JENKINS  &  ASSOC.  Columbia  SC 
FORD.  BACON  &  DAVIS,  INC  New  York  (Library! 

FRANCE  Dr  Dutcrtre,  Boulogne 

GLIDDEN  CO  STRONGSVILLE.  OH  (RSCH  LIB) 

HUGHES  AIRCRAFT  Culver  City  CA  (Tech.  Doc  Cir) 

ITALY  M.  Catront.  Milan 

KAMAN  SCIENCES  CORP  T  Cook  Colorado  Springs.  CO 
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LOCKHEED  MISSILES  &  SPACE  CO.  INC.  Sunnyvale.  CA  (K  L.  Krug) 

MCDONNEL  AIRCRAFT  CO.  Dept  501  (R.H.  Feyman),  St  Louis  MO;  Navy  Tech  Rep  V  Louis,  MO;  R 
Carson.  St.  Louis,  MO 

NEWPORT  NEWS  SHIPBLDC  &  DRYDOCK  CO  Newport  News  VA  (Tech.  Lib  ) 

NORWAY  DET  NORSKE  VERITAS  (Library),  Oslo;  DET  NORSKE  VERITAS  (Roren)  Oslo;  I.  Toss,  Oslo; 

J.  Creed,  Ski;  Norwegian  Tech  Univ  (Brandtzaeg).  Trondheim 
PORTLAND  CEMENT  ASSOC  Skokie  IL  (Rich  &  Dev  Lab,  Lib  ) 

RAND  CORP.  Santa  Monica  CA  (A.  Laupa) 

RAYMOND  INTERNATIONAL  INC  E  Colle  Soil  Tech  Dept.  Pennsauken,  NJ 

SANDIA  LABORATORIES  Library  Div.,  Livermore  CA;  Seabed  Progress  Div  453b  (D.  Talbert)  Albuquerque 
NM 

SCHUPACK  ASSOC  SO.  NORWALK,  CT  (SCHUPACK) 

SHELL  OIL  CO.  HOUSTON,  TX  (MARSHALL) 

SWEDEN  Cement  &  Concrete  Research  Inst.,  Stockholm;  GeoTech  Inst 
SWITZERLAND  T  Schneider.  Zurich 

THE  NETHERLANDS  Director  Tech  Rich  Print  Maurits  Lab,  Rijswijk 
TRW  SYSTEMS  REDONDO  BEACH,  CA  (DAI) 

UNITED  KINGDOM  Atomic  Energy  Auth  (F  R  Farmer)  London;  Cement  &  Concrete  Assoc  Wexham  Springs. 
Slough  Bucks;  D.  Lee.  London;  D.  New,  G.  Maunsell  &  Partners,  London;  R.  Browne,  Southall. 

Middlesex;  Taylor,  Woodrow  Constr  (0I4P),  Southall.  Middlesex 
WESTINGHOUSE  ELECTRIC  CORP.  Annapolis  MD  (Oceanic  Div  Lib,  Bryan);  Library.  Pittsburgh  PA 
WESTINTRUCORP  Egerton,  Oxnard,  CA 

W1SS,  JANNEY,  ELSTNER.  &  ASSOC  Northbrook,  IL  (D.W.  Pfeifer) 

WOODWARDCLYDE  CONSULTANTS  (A.  Harrigan)  San  Francisco 

BRAHTZ  La  Jolla,  CA 

BROWN,  ROBERT  University,  AL 

BRYANT  ROSE  Johnson  Div  UOP,  Glendora  CA 

BULLOCK  La  Canada 

F  HEUZE  Alamo,  CA 

CAPT  MURPHY  Sunnyvale.  CA 

R  F  BESIER  Old  Saybrook  CT 

T  W.  MERMEL  Washington  DC 
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